BACKGROUND OF THE INVENTION 
Field of the Invention 

The present invention relates to computerized methods and apparatus for 
creating particle packs and, more specifically, to computerized methods and 
5 apparatus for creating particle packs with enhanced speed and processing 
efficiency. 

Descri ption of the Related Art 

The present invention relates to the creation of particle packs that can be 
used to model or simulate particle containing or particle filled systems. There are 
J;10 many varieties of particle containing systems. Specific yet merely illustrative 
ni examples are found in the fields of asphalts, concretes, ceramics, soils, chemical 
Q physics of amorphous materials, slurry mechanics, and solid propellants or gas 
L generators. An important class of such systems comprises polymeric materials with 

particles embedded or embodied in the polymer body or binder, 
p 15 It is highly desirable in many circumstances to have a capability to model or 

simulate the properties, performance, etc. of such particle containing systems. In 
designing such particle filled systems, for example, substantial trial and error can 
be avoided if one is able to use a model or simulation to assess the implications of 
making changes in the components, their configuration, the conditions, etc. 
20 This is particularly true in microstructural analyses, e.g., wherein the 

analysis focuses on the material at a highly magnified level, such as at the 
molecular or grain stages. A specific example of microstructural analysis would 
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involve the modeling or simulation of a particle pack for a propellant, which 
typically would include fuel and oxidizer particles in a polymeric binder. 
When modeling or simulating a particle containing system, it has become 
customary to construct an analytical replica or representation of the particle filled 
5 material, which is generally referred to as a particle pack. The term "particle pack" 
has come to be used in the field to refer generally to a collection of particles in a 
defined volume. 

The general object in creating particle packs is to model or simulate the 
manner in which a plurality of particles, usually a large number of such particles, 
yjlO will fill or otherwise occupy a particular space. The nature of the particles may 
fi vary, in some cases considerably, depending on the particular application. These 
r 3 variations may include particle size variations, as well as other characteristics, 
s The complexity of models used to simulate real world particle containing 

; IJ systems increases dramatically if more than a limited number of variables are 
J=- i 15 considered. The numbers of interactions, combinations and permutations brought 
to bear by anything other than relatively simple systems causes the processing 
requirements to become unwieldy even for the most advanced computer systems. 
This complexity and processing requirement limitation is even more pronounced for 
systems that employ large numbers of particles. 
20 These limitations have caused workers in the field to invoke various 

simplifying assumptions to render the problem more workable. One such 
assumption involves the shape of the particle. The complexity and related 
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processing demands can be substantially reduced, for example, if one assumes that 
the particles are spherical. 

The general approach to simulating particle containing systems in the past 
has involved a more rigorous one in which a particle pack is constructed by placing 
5 the individual particles into a three dimensional volume, i.e., into a "full field." 
This approach in theory offers promise for accurate simulation of the real particle 
containing systems. The large number of particles required to construct a 
statistically meaningful or useful particle pack, however, usually renders this 
approach impractical. The processing power required to construct such a particle 

10 pack is prohibitive for all but extremely small numbers of particles, and only 
relatively small variations in particle size. 

The simplifying assumptions necessary to reduce processing demands down 
to a manageable level can prevent the resulting particle pack from accurately 
portraying the actual packs that are intended to be simulated. 

15 An approach that has been used to make the processing demands more 

manageable and yet the accuracy and reliability of the simulation to be good 
involves a reduced dimension approach. An example of this approach is described 
in L Lee Davis and Roger G. Carter, "Random Particle Packing by Reduced 
Dimension Algorithms," J/. AppL Phys., 67(2), 15 January 1990 pp. 1022-1029 

20 (hereinafter "Davis and Carter"). The reduced dimension approach as described 

there involves simplifying the pack construction problem by employing the principle 
that an arbitrary straight line drawn through a uniformly random three- 
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dimensional particle pack (regardless of whether the particles are spheres) will, on 
the average, lie inside particles the same fraction of its length as the packing 
fraction. As the length of the line goes to infinity, the fraction of its length residing 
inside the particles of a random pack goes to the packing fraction exactly. The line 
5 is referred to in Davis and Carter as a "central string." This approach has proven 
quite useful in enabling more complicated particle packs to be constructed. It can 
accommodate a substantially wider range of particle sizes. 

Notwithstanding the usefulness of the reduced dimension approach to 
constructing particle packs, its direct application has been limited to a certain 

10 extent in that its accuracy can be compromised by considering only particles 

intersected by the central string. In real three-dimensional packs, and assuming 
the particles are spherical, each sphere initially intersecting the central sphere, as 
it slides down the central string during placement, comes to rest more often on a 
"nearest-neighbor" sphere than on a lower central string particle. The sphere then 

15 would roll away from the central string, thereby spreading the spheres over a 
greater distance along the central string. 

A technique for addressing this limitation, also presented in the Davis and 
Carter paper noted above, involves relaxing the requirement that spheres must 
intersect the central. Under this modified approach, particles are permitted to 

20 come to rest within a defined range of the central string, after assuming that their 
path is modified to encounter and avoid previously placed spheres. This is referred 
to generally as "perturbation," or perturbed central string packing. 



This modified reduced dimension approach also has been limited, however, in 
that, when perturbation approaches are used, processing demands once again 
increase, in some cases substantially. This can give rise to a fuller and more 
complicated particle field that once again begins to approach the complexity and 
5 processing intensity of a full field. 
Objects of the Invention 

Accordingly, an object of the present invention is to provide apparatus and 
methods for creating a particle pack that simulate the actual particle distribution 
within a space with a reasonable degree of accuracy. 
10 Another object of the invention is to provide apparatus and methods for 

creating a particle pack that involve reduced processing demands and increased 
processing efficiency relatively to full field simulation. 

Additional objects and advantages of the invention will be set forth in the 
description that follows, and in part will be apparent from the description, or may 
15 be learned by practice of the invention. The objects and advantages of the invention 
may be realized and obtained by means of the instrumentalities and combinations 
pointed out in the appended claims. 

SUMMARY OF THE INVENTION 
To achieve the foregoing objects, and in accordance with the purposes of the 
20 invention as embodied and broadly described in this document, a machine- 
implemented method is provided for placing a plurality of particles, each having a 
characteristic dimension, to create a particle pack, the plurality of particles 



comprising N categories of the particles, the characteristic dimension of the 
particles of a given category being different from the characteristic dimensions of 
the particles of other ones of the categories, the characteristic dimension of the 
particles increasing as the category N increases. The method comprises a)defining 
5 a central string, a space disposed about the central string, and N concentric 
subspaces disposed about the central string and within the space, each of the N 
subspaces corresponding to one of the N particle categories; b) selecting a particle 
from the plurality of particles; c) placing the selected particle in the corresponding 
subspace so that the selected particle becomes a placed particle at a particle 

10 location unique to that placed particle and is in non-overlapping relation with other 
placed particles, the selected particle placement according to one aspect of the 
invention includes defining a catch net representative of buoyancy of a portion of 
the placed particles and positioning the catch net within the space based upon the 
placement of the portion of the placed particles. The selected particle placement 

15 according to another aspect of the invention includes defining a water level 

representative of a level of a portion of the placed particles that are smaller than 
the selected particle and represent a surface of the smaller placed particles, and 
positioning the water level within the space based upon the smaller particle surface. 
The selected particle being placed in non-overlapping relation with respect to the 

20 catch net and the water level; and d) repeating the particle selection (b) and 

placement (c) until each of the particles of the plurality of particles has become one 
of the placed particles. 



In accordance with preferred but optional implementations of the 
invention, the following may apply. The particles comprise spheres and the 
characteristic dimension of each of the particles comprises a radius. The central 
string comprises a line within the space, preferably the z axis of a rectilinear 
5 coordinate system imposed within the space. Each of the subspaces has a 
characteristic dimension that corresponds to the characteristic dimension of the 
particles of the corresponding subspace. Each of the subspaces may comprises a 
cylinder. The subspaces comprise concentric cylinders disposed about the 
central string, each of the cylinders defining one of the subspaces and having a 

10 radius that corresponds to the radius of the particles in the category corresponding 
to that cylinder. The particle selection may comprise selecting the particles 
randomly or selecting the selected particle from a distribution of the particles. 

In accordance with another aspect of the invention, an apparatus is provided 
for placing a plurality of particles, each particle having a characteristic dimension, 

15 to create a particle pack, the plurality of particles comprising N categories of the 
particles, the characteristic dimension of the particles of a given category being 
different from the characteristic dimensions of the particles of other ones of the 
categories, the characteristic dimension of the particles increasing as the category N 
increases. The apparatus comprises a) an input device for inputting particle 

20 selection information; b) a storage device operatively coupled to the input device for 
storing the particle selection information; c) a processor for selecting a particle from 
the plurality of particles, for placing the selected particle at a particle location 



unique to the selected particle in the corresponding subspace so that the selected 
particle becomes a placed particle at a particle location unique to that placed 
particle and is in non-overlapping relation with other placed particles, and for 
establishing a catch net representative of buoyancy of a portion of the placed 
5 particles and positioning the catch net within the space based upon the placement 
of the portion of the placed particles, the processor placing the selected particle 
being placed in non-overlapping relation with respect to the catch net. 

In accordance with another aspect of the invention, a machine -readable 
medium is provided for use in placing a plurality of particles, each particle having a 

10 characteristic dimension, to create a particle pack, the plurality of particles 
comprising N categories of the particles, the characteristic dimension of the 
particles of a given category being different from the characteristic dimensions of 
the particles of other ones of the categories, the characteristic dimension of the 
particles increasing as the category N increase. The machine-readable medium 

15 comprises a) machine executable instructions for defining a central string, a space 
disposed about the central string, and N concentric subspaces disposed about the 
central string and within the space, each of the N subspaces corresponding to one of 
the N particle categories; b)machine executable instructions for selecting a particle 
from the plurality of particles; c) machine executable instructions for placing the 

20 selected particle at a particle location unique to the selected particle in the 

corresponding subspace so that the selected particle becomes a placed particle at a 
particle location unique to that placed particle and is in non-overlapping relation 
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with other placed particles, the selected particle placement instructions including 
instructions for defining a catch net representative of buoyancy of a portion of the 
placed particles and positioning the catch net within the space based upon the 
placement of the portion of the placed particles, the selected particle being placed in 
5 non-overlapping relation with respect to the catch net; and d) machine executable 
instructions for repeating the particle selection (b) and placement (c) instructions 
until each of the particles of the plurality of particles has become one of the placed 
particles. 

In accordance with another aspect of the invention, an apparatus is provided 
10 for placing a plurality of particles, each particle having a characteristic dimension, 
to create a particle pack, the plurality of particles comprising N categories of the 
particles, the characteristic dimension of the particles of a given category being 
different from the characteristic dimensions of the particles of other ones of the 
categories, the characteristic dimension of the particles increasing as the category N 
15 increases, the apparatus comprises a) an input device for inputting particle 

selection information; b) a storage device operatively coupled to the input device for 
storing the particle selection information; c) a processor for selecting a particle from 
the plurality of particles, for placing the selected particle in the corresponding 
subspace so that the selected particle becomes a placed particle at a particle 
20 location unique to that placed particle and is in non-overlapping relation with other 
placed particles, for establishing a water level representative of a level of a portion 
of the placed particles that are smaller than the selected particle and represent a 
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surface of the smaller placed particles, and for positioning the water level within 
the space based upon the smaller particle surface, the selected particle being placed 
in non-overlapping relation with respect to the water level. 

In accordance with still another aspect of the invention, a machine-readable 
5 medium is provided for use in placing a plurality of particles, each particle having a 
characteristic dimension, to create a particle pack, the plurality of particles 
comprising N categories of the particles, the characteristic dimension of the 
particles of a given category being different from the characteristic dimensions of 
the particles of other ones of the categories, the characteristic dimension of the 

10 particles increasing as the category N increases. The machine-readable medium 
comprises a) machine executable instructions for defining a central string, a space 
disposed about the central string, and N concentric subspaces disposed about the 
central string and within the space, each of the N subspaces corresponding to one of 
the N particle categories; b) machine executable instructions for selecting a particle 

15 from the plurality of particles; c) machine executable instructions for placing the 
selected particle at a particle location unique to the selected particle in the 
corresponding subspace so that the selected particle becomes a placed particle at a 
particle location unique to that placed particle and is in non-overlapping relation 
with other placed particles, the selected particle placement instructions including 

20 instructions for defining a water level representative of a level of a portion of the 
placed particles that are smaller than the selected particle and represent a surface 
of the smaller placed particles, and positioning the water level within the space 
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based upon the smaller particle surface, the selected particle being placed in non- 
overlapping relation with respect to the water level; and d) machine executable 
instructions for repeating the particle selection (b) and placement (c) instructions 
until each of the particles of the plurality of particles has become one of the placed 
5 particles. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The accompanying drawings, which are incorporated in and constitute a part 
of the specification, illustrates a presently preferred embodiments and methods of 
the invention and, together with the general description given above and the 
10 detailed description of the preferred embodiments and methods given below, serve 
to explain the principles of the invention. 

Fig. 1 is a pictorial diagram of a computer system which comprises a 
presently preferred embodiment of the apparatus according to the invention, and 
upon which the presently preferred version of the inventive method may be carried 
15 out. 

Fig. 2 is an illustrative example of a particle pack construction according to 
the preferred implementation, wherein spherical particles are placed along the 
central string. 

Fig. 3 shows another illustrative example of a particle pack construction 
20 according to the preferred implementation, wherein spherical particles are placed 
along the central string, but wherein the particles extend more several particle 
diameters out from the central string. 
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Fig. 4 shows a cylindrical space and cylindrical subspaces within it according 
to the preferred embodiment and the preferred method of the invention. 

Fig. 5 shows an illustrative particle pack constructed using the preferred 
implementation. 

5 DETAILED DESCRIPTION OF THE 

PREFERRED EMBODIMENTS AND METHODS 

Reference will now be made in detail to the presently preferred embodiments 
and methods of the invention as illustrated in the accompanying drawings, in which 
like reference characters designate like or corresponding parts throughout the 

10 drawings. It should be noted, however, that the invention in its broader aspects is 
not limited to the specific details, representative devices and methods, and 
illustrative examples shown and described in this section in connection with the 
preferred embodiments and methods. The invention according to its various aspects 
is particularly pointed out and distinctly claimed in the attached claims read in 

15 view of this specification, and appropriate equivalents. 

In accordance with one aspect of the invention, an apparatus is provided for 
creating a particle pack. In accordance with another aspect of the invention, a 
machine-implemented method is provided for placing a plurality of particles to 
create a particle pack. 

20 A particle as the term is used herein is used according to its common but 

broad meaning in the field. It includes any type of particle that might be found in 
particle rigid body of arbitrary shape containing materials. In the preferred 
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implementation, particles comprise rigid or nearly rigid objects of arbitrary- 
It is useful in analyzing particle packs to be able to segregate the particles 
into categories. A group of substantially identical particles, or a group of particles 
that for purposes of analysis may be deemed identical, and which group is in some 
5 way distinguishable from other groups of particles is called a "mode," If there is 
only one type of particle comprising the pack and those particles are identical to 
each other, or may be assumed to be identical, the pack is said to be "monomodal." 
If there are two types of particles, the pack is "bimodal," and if multiple types of 
particles are present, it is "multimodal." 

10 A mode represents a principal size or type of particle. In practice, particles 

even within a mode often are not exactly identical, but instead have slight 
variations. In some instances, for example, each mode can be divided into a size 
distribution about some mean. The distribution may comprise a number of discrete 
sizes about that mean. The sets of various particles of a given size or character 

15 within a particular mode are called "submodes." In the presently preferred 

implementation, if there is no size distribution attached to a particular mode, the 
submode arrays simply contain the mode data as their first and only entry. 
Because both modes and submodes theoretically merely involve segregating 
particles into categories based upon one or more distinguishing characters, modes 

20 and submodes are referred to generally herein as "categories," 

The apparatus according to the invention preferably comprises a computer 
system 100, such as a commercially available personal computer, small business 
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computer, engineering work station, mainframe, microprocessor, or virtually any 
other machines with similar processing power, appropriately configured and 
programmed with computer software or "code" as outlined and described herein for 
performing the tasks identified herein. It includes a processor 102, storage in the 
5 form of RAM 104, a hard drive 106, a compact disk (CD) drive 108, and a drive for 
diskettes 110, all contained within a cabinet 112. System 100 also includes a 
monitor 116, a keyboard 118, and a pointing device 120 such as a mouse, track ball 
or the like. System 100 also may be part of a network, in which case it may include 
a network connection 122 and network bus 124. The method according to the 
10 present invention, in its presently preferred implementations, may be carried out on 
computer system 100, using the same hardware and software, as described herein 
with respect to the preferred embodiment of the apparatus according to the 
invention. 

"Machine" as the term is used herein refers broadly to a computer, processor, 
15 microprocessor, a network of such devices, or other apparatus capable of performing 
processing as described herein. Such machines typically and preferably will 
comprise a computer, such as computer system 100. 
OVERVIEW OF THE PREFERRED IMPLEMENTATION 

The method and apparatus according to the presently preferred yet merely 
20 illustrative versions of the invention will be described in this document as involving 
the implementation and execution of a computer program that can be run on a 
computer system or machine as described above. The preferred versions of the 
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apparatus and methods according to the invention are referred to collectively in this 
document as the "preferred implementation." 

The architecture of the preferred implementation is outlined below. In its 
presently preferred form, it is written in double precision C++ or Fortran 90. This 
5 is not, however, limiting. The code may, for example, be written in FORTRAN, or 
other suitable programming languages. An example of such a computer system is 
shown in Fig. 1. The code in its preferred version builds a single particle pack and 
outputs user-specified pack microgeometry and statistics. It is modular and can be 
disassembled into its component subroutines for use elsewhere. 

10 In accordance with the apparatus aspect of the invention, an input device is 

provided for inputting particle selection information, and a storage device 
operatively coupled to the input device for storing the particle selection information. 
The information used by the preferred implementation will be described 
throughout this document. The input device may comprise any one or any 

15 combination of several input means, including a network download, the CD reader, 
a diskette reader, the keyboard, the tracking device, a tape drive, or any other 
means used for inputting data into a machine capable of implementing these 
aspects of the invention. The storage device also may comprise any suitable storage 
means, and may include RAM, the hard drive, an optical storage device, storage 

20 devices on a machine operably coupled to or accessible by computer 100, or other 
suitable storage devices readable by a computer directly or with translation. 
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The main program of the preferred implementation comprises seven 



subroutine calls. The names and functions of these subroutines are as follows: 



SUBROUTINE 
NAME 


SUBROUTINE FUNCTION 


INITERROR 


Initializes arrays containing the number of errors, 
warnings, and notations accumulated during program 
execution. It also initializes arrays containing 
information on how these data are to be displayed. 


INNAME 


Looks for and reads an input file named INNAME. INP 
which contains the names of files containing the input 
data needed by the preferred implementation to build 
the user-specified pack and output the requested pack 
information. If the file INNAME.INP does not exist, 
default file names are used to find instructions on pack 
building. Default values are used if user-specified 
instructions are lacking. If essential instructions are 
lacking, pack building is aborted and error messages are 
issued. 


SETHISTO GRAM 


Reads the information for the types, quantities, and 
sizes of particles to be used in constructing the pack. 


SETCONTROL 


Reads flags and parameters controlling the construction 
of the pack. Default values are used if files containing 
these instructions cannot be found. 


SETOUTPUT 


Reads flags and parameters controlling the type of pack 
data to be output. Default values are again used if files 
containing these instructions cannot be found. 


PPACK 


Constructs the particle pack according to instructions 
read by the previous subroutines. Each time it is called 
to build a new pack, it initializes all variables and arrays 
so there will be no interfering data left over from a 
previous call. 


OUTSTAT 


Outputs the user-specified information on pack 
microgeometry and statistics. 
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In summary, the first subroutine (INITERROR) initializes error, warning, 
and notation arrays and parameters. The next four subroutines (INNAME, 
SETHISTOGRAM, SETCONTROL, SETOUTPUT) and the seventh subroutine 
5 (OUTSTAT) control the reading of building instructions and the output of data. 

The sixth routine (PPACK) constructs the particle pack. This modularity facilitates 
incorporation of the code or portions of it into larger calling codes, such as codes 
using numerical methods or finite element analysis. 

Because in the preferred implementation a substantial amount of 
10 information may be used to govern particle pack construction, the manner in which 
the data used by the software code is categorized in the context of the preferred 
implementation will now be described. These data categories are placed into 
common blocks and passed from one code section to another as needed. In the 
presently preferred implementation, all data about the pack along with building 
15 instructions are placed into one of the following categories: 



DATA BLOCK 


DESCRIPTION 


HISTO 


Governs the types and quantities of particles to be used 
in constructing the pack, i.e., the particle histogram. 


CONTROL 


Contains all parameters and flags controlling how the 
pack will be built and stored. 


CURINFO 


Contains information on the status of the particle 
currently being placed, such as its position and which 
other particles it is currently touching. 
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PLACED 


Contains information on all particles placed so far in 
the pack including the position, size, and type of each 
particle. 


SETOUT 


Contains the flags and parameters governing what pack 
statistics are to be generated and how the pack 
microgeometry and its statistics are to be output. 


ERRORS 


Contains the names of errors, warnings, and notations 
that can occur in program input, execution, and output. 
It is used to monitor the number of times each error, 
warning, or notation occurs and records such 
information in a file and/or to the screen, according to 
the dictates of the user. 



THE SUBROUTINE INITERROR 

In the presently preferred implementations, there are three categories of 
information, called "advisors," used to monitor the health of the pack, i.e., to 
5 monitor whether any of the pack construction rules, such as the prohibitions 

against particle overlap, and general requirements for particle distribution in the 
pack (e.g., homogeneity at the macro level of the random or other specified 
statistical distributions at the microphone mount level) have been violated. They 
are contained in the program's error arrays and parameters. The three categories 
10 of advisors are errors, warnings, and notations. 

Error advisors are of two kinds, fatal and nonfatal. A fatal error will stop 
program execution because there is not enough information available to build the 
pack. Nonfatal errors will not stop the program execution but will use default 
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values in place of erroneous user-data. Therefore, a pack will be built but may not 
be the kind of pack the user intended. 

Warning advisors are notifications of unusual situations or potential 
problems. If the warning deals with user input, the user should correct the input 
5 problem because the warning sometimes precipitates the use of defaults that may 
not have been intended. 

The third category is notation advisors. There are many data that may be of 
interest during pack building. Notations report on these items. They are not errors 
or warnings but simply notifications of items of interest as pack building proceeds 
10 or after it has finished. 

All advisors in this implementation have arrays that keep track of how many 
times each has been called: NumErrors(0 ? NumWarnings(i), and NumNotes(i), for 
the i-th error, warning, or notation, respectively. The totals for each category of 
advisor are also accumulated in the parameters NTotErrors, NTotWarnings, 
15 NtotNotes. 

The display of advisors in this preferred implementation is controlled by 
arrays containing integers limiting the number of times an advisor is to be 
displayed to the computer screen and written to a disk file created for that purpose. 
The integer arrays LimScErrors(i) and LimDiErrors(i) dictate how many times the 
20 i-th error is to be written to the screen and the disk file, respectively. Similarly, 
LimScWarnings(i) and LimDiWarnings® are provided for warnings and 
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LimScNotes(i) and LimDiNotes® for notations. If any of these array elements is 
less than or equal to zero, the corresponding advisor is not written. If the user 
wants all advisors written every time they are encountered, the limit arrays need 
simply be set to a very large integer (in this implementation not to exceed 
5 2147483646, which is one less than the maximum integer size). The total number 
of times any advisor has been encountered is monitored in NTotErrors, 
NTotWarnings, and NTotNotes whether or not it is displayed. 

All data on all advisors according to this preferred implementation are kept 
in and monitored by the subroutine ErrWarnNote(ierr, iwarn, inote, iunitnumber), 

10 where ierr is the integer specifying the error, and similarly for iwarn and inote. The 
integer iunitnumber is the disk file unit number to which advisors are written as 
requested. It is assigned the number 13 in this implementation. 
THE SUBROUTINE INNAME 

In the presently preferred implementation, the data required by the code are 

15 read from several data files. The INNAME subroutine reads the names of the data 
files in which the data resides from a file called INNAME.INP, The files whose 
names are specified in INNAME.INP are as follows: 



FILE 


DESCRIPTION 


DEFAULT NAME 


Histogram 


Contains data on sizes, quantities, and types of 
particles along with material information and 
the sizes of the cylinders into which they are to 
be dropped. 


HISTJNP 


Drop File 


Contains data on where and how particles are to 
be dropped onto the pack. 


DROPDATA.INP 
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Preplace File 


Containing the positions and types of particles to 
be preplaced in the pack, if any. 


PREPLACEJNP 


Input Control 
File 


Contains the input flags and parameters that 
control how the pack is to be built. 


INFLAG.INP 


Output Control 
File 


Contains flags and parameters controlling what 
output data is to be generated and how it is to be 
displayed. 


SETOUT.INP 



In accordance with the preferred implementation, if the file INNAME.INP 
does not exist, then default names are used for each of the five files containing input 
data. The default names for this implementation are identified in the third column 
5 of the table above. 

In the preferred implementation, the histogram file is mandatory but the 
other files are not. Packs can be built if all other files are missing by using defaults. 
Of course, the various options and data production capabilities of the presently 
preferred implementation may not be exerciseable without these other files. 
10 THE SUBROUTINE SETHISTOGRAM 

The SETHISTOGRAM subroutine according to this preferred implementation 
reads the particle histogram from the histogram file and sets necessary parameters 
in the program associated with that file. The particle histogram comprises a set of 
particle modes and submodes, in this implementation in the form of arrays. If there 
15 is no size distribution attached to a particular mode, the submode arrays simply 
contain the mode data as their first and only entry. 
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First, subroutine SETHISTOGRAM checks to see if the histogram file exists. 
If it does not, the program is terminated with a fatal error message. If the file 
exists, the data is reads as follows. 

The number on the first line of the histogram file is the "stop build" criterion 
5 in the form of an integer called STOPBUILD. If STOPBUILD is positive, then that 
is the total number of particles to be placed in the pack. If STOPBUILD is less than 
or equal to zero, then pack building is to be stopped when the center of the highest 
placed particle reaches a certain altitude or when a pre-specified maximum number 
of particles is placed, whichever comes first. If STOPBUILD is less than or equal to 
10 zero, the second line is the real number STOPALTITUDE specifying the altitude at 
which pack building is to cease. 

The particle histogram is processed next in this implementation. Each set of 
lines contains the data for one mode. The modes do not have to be in any particular 
order. The first line of the set contains the mode information. If there is a size 
15 distribution associated with this mode, a set of lines follows that specifies the 
submodes. In this implementation there are six numbers on each mode line. 

1. The first number is the mean particle size for the mode. It is read in as a 
particle diameter but is internally converted to particle radius. 

2. The second number is the mass quantity for the mode. Because the total 
20 mass quantities are normalized later, relative mass quantities are adequate. 

3. The third number is the mass density for the mode. The units of particle size 
and density preferably are consistent. In this implementation this is a 
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requirement. For example, one should not input a mode diameter in microns and 
its density in grams/cm 3 . 

4. The fourth number is an integer specifying a material number for the mode. 
Material numbers are a way of specifying some property unique to a mode or set of 

5 modes. For example, a material number may specify a particular material 
composition or it may specify a set of modes that have some other feature in 
common. 

5. The fifth number is the radius of the cylinder into which the particles of this 
mode will be dropped. The centers of the particles are constrained to stay within 

tfO this cylinder wall. If there is a distribution associated with this mode, the cylinder 
%l s radius of the mode is ignored and the cylinder radii for each submode are used, 
fn 6. The last number is a logical parameter RDIST specifying whether there is a 

P distribution of submodes associated with this mode. If the parameter is FALSE, 

there is no distribution and the first submode array elements are filled with 
If information for that mode. If the parameter is TRUE, the program next reads the 
M submode information in the following form. 

The submode data for this implementation comprises a set of lines, each one 

of which includes the following three pieces of information: 

1. The submode size. As with the mode sizes, it is read in as a diameter but 
20 internally converted to a radius. 
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2. The fractional mass of the submode. This number is the fraction of the 
mode's mass contained within this submode. The sum of a mode's fractional masses 
should add to unity. The program normalizes them if they do not. 

3. The submode cylinder radius. 

In many parts of the code as presently preferred and implemented, it is 
convenient to have the submodes ordered in increasing particle size. A global size 
array is created that contains all submodes in increasing particle size. The global 
submode arrays comprise an array of submode radii (QSIZE), submode mass 
(QMASS), submode mass density (DENSITY), submode material number 
(MODEMAT), and submode cylinder radii (SUBCYLRAD). 

The end of a submode is signaled by a set of zeroes everywhere that a 
nonzero number was expected, and the end of a mode is likewise signaled when a 
set of zeroes is read. 

After reading all mode and submode information, the SETHISTOGRAM 
subroutine reads the information in the drop file if it exists. If it does not exist, the 
program drops the particles randomly anywhere within the cylinder corresponding 
to that mode or submode. If it does exist, then it contains the following data. 

The first line has two integers. The first integer is a drop option, either 1, 2, 
or 3. The second integer is the number of particles that are to be dropped according 
to that option. 

Option 1 means that first group of particles are to be dropped randomly 
anywhere within their cylinder. 
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Option 2 means that the particle group is to be dropped within a circle of 
relative radius a r where the drop radius is some fraction of the cylinder radius (0 < 
ctr < 1) of the particle that will be selected. The center of this circular region is at 
relative radial position p r where p r is some fraction of the particle's cylinder radius 
5 when it is chosen (0 < p r < 1). Because dropping particles outside the cylinder wall 
is forbidden, it should hold that p r + a r < 1. These data are stored in the pack 
position arrays. 

Option 3 means that each particle in that group is to be dropped from a 
specified position and is of a specified mode. The x position of the drop site is stored 

CI 10 in the x position array, the y position in the y position array, and the mode number 
in the mode number pack array. If an (x,y) position lies outside the cylinder 

?^ corresponding to that mode, a warning message is issued and the particle is 

e randomly dropped anywhere within the cylinder. 

~* The end of this input file is signaled by a value of zero for the drop option. 

JL;fL5 Another function of the SETHISTOGRAM subroutine in this implementation 

is to create an index array that orders the global cylinder radii in ascending order 
with an integer array called CylSizeOrder so that cylrad(CylSizeOrder(l)) is the 
smallest cylinder radius, cylrad(CylSizeOrder(2)) is the second smallest, and so 
forth. 

20 THE SUBROUTINE SETCONTROL 
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The currently preferred implementation includes a variety of options for pack 
construction. Some of those options are read in a file that can have any name as 
specified in Subroutine INNAME but whose default value is INFLAG.INP. This 
subroutine SETCONTROL reads those control parameters. 
5 If a data file containing input control parameters exists, then the control 

parameters are read from it in this implementation. If it does not exist, a default is 
set for each parameter in this routine. The parameters used in this implementation 
are described in detail below but in summary include the following: 

ROUGHNET = A flag that specifies whether the k-th mode catch net 

10 (described below) is ratcheted up after each particle placement so that no two 

particles are at exactly the same altitude. It is mostly used at the beginning of pack 
building so the spheres at the bottom of the pack do not lie exactly on the bottom of 
the cylinder. The default is TRUE. 

POOL = A flag that specifies whether the number of particles in the pack are 

15 selected by a random number generator, weighted by their frequency of occurrence 
in an infinite pack, or whether they are selected from a pool whose mode 
proportions are, to the nearest integer, exactly the proportions of an infinite pack. 
If the POOL flag is TRUE, then random chance variations in modes that contribute 
few particles to the pack are avoided. The default is TRUE. 

20 FINDLOW = A flag that orders find the lowest exposed north pole of placed 

particles within its drop cylinder and to drop immediately in its neighborhood. It is 
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overridden by any instructions that may appear in the drop file. The default is 
TRUE. 

NDROP = A positive integer that specifies the number of times a trial drop of 
each particle is to be made. The lowest final resting place is selected from among 
the trial drops. Higher values for NDROP means the pack will be more dense. The 
default is one. 

ISEED = An integer specifying the random number seed. A non-positive 
integer will trigger to seed the random number generator from the system clock. 

TUNNEL = A flag that gives permission to check for cavities underneath 
larger spheres when a smaller sphere is dropped. The small sphere "tunnels" 
through the large one and, if there are no overlaps in the southern or underneath 
side, is dropped again from that underneath position whereupon it searches for a 
final resting place in that cavity. The tunneling procedure seeks to mimic tamping 
or shaking a pack to fill in such cavities, thereby making the pack more dense. 
THE SUBROUTINE SETOUTPUT 

A particle pack is typically created so that the system or method user can 
obtain specific information about the actual particle pack the simulation is intended 
to simulate. The information required depends on the application. The subroutine 
SETOUTPUT in this implementation reads the flags and parameters that instruct 
other parts of the preferred implementation to create specific pack information. In 
this implementation, for example, if the instruction file for setting output flags and 
parameters does not exist, the code outputs a limited amount of information into a 
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file called PACKDATA.OUT. Otherwise, this basic pack information is output to a 
file name specified in Subroutine INNAME. With no specifications, the following 
data is sent to the output file: 

• an echo of the particle histogram; 
5 • the packing fraction; 

• the maximum and minimum positions of particles for each submode as well as 
overall maxima and minima; and 

• the number of particles placed for each submode as well as total particles placed. 

If the pack output file exists when using the preferred implementation, it 
10 preferably contains the following information. The first line contains a flag 

specifying whether the particle pack is to be written to a disk file. The default is 
FALSE. If the flag is TRUE, then the next set of lines specify how the pack file is to 
be written. On the first subsequent line, the format of the pack (ASCII or binary) is 
written. On the second subsequent line, a flag is read that determines whether the 
15 material number and radius are written in the pack file in addition to the three 
coordinates and the mode number. On the third subsequent line, the output 
geometry (Cartesian or cylindrical coordinates) is written. The name of the file to 
which the pack is to be written is on the fourth line. 

The next line contains a flag specifying whether radial distribution functions 
20 gij(r) are to be generated. If it is TRUE, then the functions that are to be generated 
are specified in a list of integer pairs i and j on the subsequent lines, where the i 
and j represent modes or submodes. Each line also contains a Ar which is the 
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increment in terms of the j-th particle radius in which the number of particle 
centers are to be counted (default = 0.1 radius) and rmax which is the extent to which 
the radial distribution function extends in terms of j-th particle radius. Hence each 
radial distribution function requested is on a line with the format 
5 i j Ar rmax . 

The list ends when zeroes are read. Each radial distribution function is written to a 
file with the name g[i]_[/].dat where the square brackets indicate the actual 
specified integers to be used in the file name, e.g., #3_14.dat for the radial 
distribution function of Mode 14 spheres about Mode 3 spheres. The default for the 

10 radial distribution flag in this implementation is FALSE. 

The next line in this implementation contains a flag specifying whether 
coordination numbers are to be generated. If it is TRUE, then the coordination 
numbers that are to be generated are read as a list of integer pairs i and j on 
subsequent lines. The list ends when zeroes are read. The coordination numbers 

15 preferably are all written to one disk file COORNUM.DAT. The default for the 
coordination number flag in this implementation is FALSE. 

Radial distribution functions and coordination numbers and their calculation 
procedures are defined in the discussion of subroutine OUTSTAT, below. 

In accordance with the apparatus aspect of the invention, the apparatus 

20 comprises means for defining a space, a plurality of subspaces within the space, and 
a central string within the space and within the subspace. In accordance with the 
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method aspect of the invention, the method comprises defining a space, defining a 
plurality of subspaces within the space, and defining a central string within the 
space and within the subspaces. 

The subspaces preferably comprise N concentric cylinders disposed about the 
5 central string and within the space. Each of the N cylinders corresponds to one of 
the N particle modes and has a cylinder radius corresponding to the particle radius 
of the particles of the corresponding particle mode. The means as implemented in 
the preferred embodiment comprises the processor. 

The space and subspaces are constructions that define the regions or volumes 

10 into which the respective particles are to be placed. The space may comprise any 
contained geometric region, but preferably comprises a cylindrical construction. 
The subspaces also may comprise any contained geometric region, but preferably 
have a shape and configuration that corresponds to the overall space. In the 
preferred versions of the apparatus and method, the space comprises a cylindrical 

15 volume, and the subspaces comprise a plurality of concentric cylinders, disposed 
about a longitudinal or z axis, wherein sequential ones of the cylinders have a 
radius that is larger than the prior cylinders. An illustrative example of such a 
concentric cylinder configuration is shown in Fig. 2. The number of cylinders in the 
preferred implementation is equal to the number of particle categories, e.g., 

20 submodes, that will be present in the particle pack. The radius of each cylinder is 
selected to correspond to the radii of the particles. The radius of the first cylinder 
corresponds to or is a function of the particle radius of the smallest mode, the radius 
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of the second cylinder is equal to the particle radius of the second smallest mode, 
and so on. The cylinder radius for a particular particle category preferably 
corresponds to or is a function of the radius or a representative radius of the 
corresponding particle for that mode or submode. The cylinder radius may be any 
5 positive number multiple of the particle radius. Preferably but optionally the 

cylinder radius is about 3 to 10 times the particle radius, and move preferably about 
5 to 10 times the particle radius. 
THE SUBROUTINE PPACK 

Further in accordance with the method aspect of the invention, the method 
10 includes selecting a particle from the plurality of particles, and placing the selected 
particle in the corresponding subspace so that the selected particle becomes a placed 
particle at a particle location unique to that placed particle and is in non- 
overlapping relation with other placed particles. 

In accordance with the preferred implementation, particles are selected 
15 randomly from the distribution of particles using a Monte Carlo method, although 
this is not necessarily limiting. In terms of the apparatus, this comprises a random 
number generator in computer 100. 

Regarding terminology, a "selected particle," also referred to herein as a 
"current particle," is a particle that has been selected for placement, but has not yet 
20 been placed (become a placed particle") in a particle location in the particle pack. A 
"placed particle" is a particle that has been placed or fixed by assigning it a unique 
particle location in the particle pack. 
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In the preferred implementation, the software code subroutine selects 
particles and constructs the particle pack according to instructions read by the 
previous subroutines. Each time it is called to build a new pack, it initializes all 
variables and arrays so there will be no interfering data left over from a previous 
5 call. Code lines of PPACK and a brief explanation of each will now be provided. 
CALL PREPAREPACK 

Subroutine PREPAREPACK performs three functions in this 
implementation. First, it preplaces particles, if any, that have a position assigned 
by the user, such as would be desirable if semi-random packs were to be built. 

10 Secondly, it arranges histogram data such that particles can be chosen by the 

random number generator with probabilities appropriate for their abundance in the 
reduced dimension pack. Third, it initializes and sets parameters that will be used 
by the rest of the subroutines in PPACK. 
DO icarus = 1, stopbuild 

15 This DO loop continues to cycle until either the user-specified number of 

particles have been placed in the pack or until a pack height has been reached. 
CALL CHOOSE 

Subroutine CHOOSE chooses the type, category, mode or submode of the 
next particle to be dropped. It drops the "current particle" according to user- 
20 specified instructions. 
DO idrop = 1, ndrop 
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The user may specify more than one drop at different drop sites for each 
particle so that the lowest final site can be selected. Multiple drops are equivalent 
to tamping the pack so the particles settle into lower sites resulting in a denser 
pack. 

5 CALL DROPSITE 

Subroutine DROPSITE selects an (X,Y) position from which the particle will 
be dropped. The site can be chosen randomly, or the user may ask the code to use a 
criterion to find potential sites where the particle is likely to settle to an especially 
low final resting place. 

10 CALL PLACE 

Subroutine PLACE finds the final resting place or "placed particle location" 
for the current particle, which is the site at which the current particle will 
ultimately be placed. The subroutine PLACE in this implementation contains all 
the instructions on how a particle drops, rolls, and encounters other objects in the 

15 pack as it moves toward its final resting place. It is where the machine spends the 
great majority of its processing time during pack construction. 
END DO 
CALL UPDATE 

Subroutine UPDATE updates all the particle pack arrays and parameters 
20 with the position and additional data associated with the newly placed particle. 
END DO 
RETURN 
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END 

In accordance with the preferred implementation, a particle pack is 
constructed by "dropping 1 ' the selected particles one at a time into the set of 
concentric cylinders. The dropping presumes a vertical orientation, which is not 
5 necessarily required. But regardless of the reference frame, the current particle is 
moved in the space until it encounters another object, as described in principle 
herein. Each particle mode is dropped somewhere within the cylinder radius of its 
corresponding cylinder. The radius of the cylinders in this preferred 
implementation is specified by the user, but preferably the small particles are 

10 dropped into a small-radius cylinder and the larger particles into a larger-radius 
cylinder or subspace. The large particles can therefore overlap the inner small- 
radius cylinders but the small particles do not have to fill the volume in the larger 
cylinders. The reduced dimension approach as described in Davis and Carter is 
used in this preferred implementation. 

15 The reduced dimensional approach can make pack building tractable from a 

processing demand perspective when the particles comprising the pack have broad 
size ranges. Particle centers in this implementation cannot pass outside the radius 
of the cylinder wall to which they belong but can pass through any smaller cylinder 
walls. Because only the innermost cylinder can have all modes or submodes 

20 intersecting it, only there does the reduced dimension pack literally approximate 
the corresponding real three-dimensional pack. The remainder of the pack can 



34 



however be simulated using this approach, for example, by assuming that the 
remainder of the space is statistically identically configured. 

In the presently preferred implementation, the selection of particles to be 
placed and the placement of the particles in the particle pack are undertaken by the 
5 processor within the machine under the execution of the PPACK subroutine. In 
this implementation, there are four subroutines called by PPACK 
(PREPAREPACK, CHOOSE, PLACE, and UPDATE). Each of these subroutines 
and the functions they perform will now be described. 
THE SUBROUTINE PREPAREPACK 
10 As noted above, the PREPAREPACK subroutine has three functions: 

preplacing particles, setting cumulative arrays so the random number generator 
can readily choose the mode of the next particle to be placed, and initializing 
various parameters in preparation for building the pack. Each of these functions 
will now be described in turn. 
15 Pre-Placing Particles 

The first function of PREPAREPACK is to pre-place any particles for which 
the user has specified a position. This option permits semi-random or fully ordered 
packs to be constructed. A filename is read in Subroutine INNAME that contains 
the global (X, Y y Z) position and the mode or submode number of each particle to be 
20 pre-placed. If the filename is not read, the default filename for preplaced particles 
is PREPLACE.INP. If a filename of preplaced particles does not exist, the program 
assumes no particles are to be preplaced and this option is not exercised. The 
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preplace file contains five numbers for each particle to be placed. In order, they are 
the X position, the Y position, the Z position, and the mode number, and the 
submode number. These are read into the surface arrays. 
Setting Arrays For Choosing Particles 
5 The second function of PREPAREPACK according to the preferred 

implementation is to make an array so that the random number generator can 
randomly choose, with the appropriate probability, the submode of the next particle 
to be placed in or on the pack. The relative fraction of particle numbers for each 
submode therefore should be calculated. The number of particles in each submode 

10 to be placed will depend on the submode mass, its mass density, its particle radius, 
and its cylinder radius, e.g., in the following manner. 

In this implementation there are two options in setting up the weighting 
factors that govern how the random number generator will choose the next particle 
to be placed. The first option is to use the same weighting factors each time a 

15 particle is chosen for placement. This is referred to herein as the "unchanged- 
weighting option." The second option, referred to herein as a "pool option," creates a 
pool that has a set of bins with the appropriate number of particles placed into it. 
Each time a particle is drawn from one of the bins, the weighting factor of drawing 
from each bin is altered slightly. Creating a pool of particles from which to draw 

20 forces the right proportion (to within nearest integers) of particles from each 

submode to be placed. In the limit of an infinite number of particles in the pack, 
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these two methods are identical. The pool option more closely resembles a very 
large pack. The unchanged-weighting option will be described first. 
Unchanged-Weighting Pool Option 
The total volume of k-th submode particles to be placed is 



4 3 



(1) 



where mk is the submode mass, pk is the submode mass density, and ak is the 
particle radius for that submode. 

The relative number of submode k particles in this reduced dimension pack is 
also proportional to the area of the confining cylinder for that submode: TiRk 2 where 
10 Rk is the radius of the cylinder. Therefore, the fraction fk of fe-th submode spheres to 
be placed in a reduced dimension pack is proportional to 

(2) 

Pk a k 

To obtain the actual fraction fk, the above expression can be normalized: 
F = t& (3) 

k=\ 

15 where N is the total number of particle submodes. Then 

fk=<pk/F. (4) 

These are the unchanged-weighting factors that determine which particle submode 

shall be chosen for the next placement. 
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These weighting factors are placed into a cumulative array so the random 



number generator can immediately choose the submode. In other words, the 



cumulative array is the sum of all weighting factors prior to it: 



Ci 



fi 



5 



(5) 



Cn=1. 

When a uniform random number r is generated on [0,1], a bisection search is made 
of the array {Ci,C2,...Cn} to find which pair of numbers the random number lies 
between. The k-th submode is chosen if Ck-i <x<Ck. Co is defined to be zero. 



The pool option is similar. As discussed before, the user can specify the 
number of particles to be placed in the pack or can specify a pack height in which 
case the program continues to add particles until the height is achieved (or until the 
array sizes are exhausted). The pool option in this implementation can only be used 
15 if the user has specified the total number of particles to be placed in the pack, say n. 
The number of submode k particles nk to be placed is 



where NINT indicates the integer nearest to its argument. Because of nearest 
integer roundoff, it is possible that the nk do not sum exactly to n. Then the nk 
20 quantities should be adjusted to compensate where the roundoff errors were the 
largest. 
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Pool Option 



n k = NlNT[nfk] 



(6) 
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It is possible, especially for packs with very broad particle size distributions, 
that some nfk < % so that submode would not be represented at all in the pack. The 
preferred implementation issues a warning when that happens, or even when a 
particular m is very small, so that the user is aware that the pack is too small to be 
5 able to sample enough different configurations to mimic a true three-dimensional 
particle pack. Changing cylinder sizes on some submodes can solve this situation, 
should it arise. 

INITIALIZING VARIOUS PARAMETERS 
The packing routine can be called many times, as when the user wishes to 
10 search for an optimum pack subject to some set of constraints. Some of the 

parameters and arrays from previous calls must not be carried into successive calls 
or they can taint the results. This portion of Subroutine PREPAREPACK 
initializes all necessary parameters, flags, and arrays for use in a new pack 
simulation. 
15 THE SUBROUTINE CHOOSE 

The CHOOSE subroutine uses the random number generator to choose the 
submode of the next particle to be added to the pack. The random number 
generator in this implementation generates a uniform distribution on [0,1]. The 
cumulative array from the previous routine PREPAREPACK is a set of numbers 
20 such that 

Co = 0 < Ci < C 2 < ... Cn-j <Cn=1. (7) 
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Because this cumulative array can be large, and because there is no 
restriction on the relative distances between successive array elements, in this 
implementation, a search by bisection is a preferred method to find two numbers 
between which a given random number lies. 
5 When a random number r is chosen in [0,1], we want to find the two array 

elements that bracket it such that 

Ck-i<r<Ck , l<k<N. (8) 

(In the unlikely event that r = 0, it will be assumed to belong to the first bin, i.e., 

where k = l.) In preferred bisection methods, a check is made to determine on 

10 which side of Cn/2 the random number lies. Then attention is concentrated on that 
side and it is bisected again, continually narrowing the range of {Ck} in which r lies 
until Eq. (8) is satisfied. 

As noted above, the particle just chosen is referred to herein as the "current 
particle" ox the "selected particle" in the discussion below. 

15 THE SUBROUTINE DROPSITE 

The DROPSITE subroutine chooses a drop site from which the current 
particle can be dropped. Actually, if the user so specifies, it chooses several drop 
sites and then keeps the site at which the particle came to the lowest possible final 
resting place, referred to herein as a "placed particle location." In a sense, this is 

20 the same as tamping the pack so that it settles into a denser configuration. The 
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user specifies how many trial drops a particle is to undergo. The default in this 
implementation is one. 

There are several options available to the user in this implementation when 
deciding where a new selected particle is to be dropped onto the pack. 
5 1. Drop it randomly anywhere within its cylinder; 

2. Drop it randomly anywhere within a user-specified circle that is contained 
within the cylinder; 

3. Drop it from a user-specified (X,Y) position; and 

4. Select the lowest exposed north pole among the placed pack spheres and drop 
10 it randomly within a short distance of that pack sphere's ( X, Y) position. If dropping 

multiple times, it is preferable to select the several lowest exposed north poles and 
drop within a short distance of each of them. 

Dropping Randomly Within The Cylinder 
When dropping a sphere anywhere within a given cylinder, two random 

15 numbers are generated, u 9 and u p . The azimuthal angle (in radians) of the drop 
position <jk relative to the central string is chosen by 
(pc - 2%u (p . (9) 
The radial position is not chosen uniformly because there is more area for large 
radial positions in a circle than for small radial position. Because the area grows as 

20 the square of the radial size, all areas will be equally accessible if we take the 
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square root of the random number generator and multiply it by the cylinder wall 
radius for the current sphere W c : 

P c =Kfe~ P - do) 

The (X, Y) drop site is found as follows: 

5 Xc- pe COS^c (11^) 

Yc = pc sin^c . (lib) 

Dropping Randomly Within A User-Specified Circle 

Suppose the user specifies that the particles are to be dropped randomly 
10 within a circle of radius R centered at (Xo, Yo). In the preferred implementation it is 
required that the circle be contained within the cylinder of the current sphere, i.e., 



^X 2 0 +Y 0 2 +R<W C . (12) 

Otherwise, the circle is an illegal specification because a particle cannot be dropped 
outside its cylinder wall. In that case, the preferred implementation will default to 
15 dropping randomly anywhere within the cylinder wall at Wc. A value of zero for R 
is permitted, meaning the particles are dropped from the position (Xo, Yo). If this 
condition is satisfied, then the drop site can be selected as follows. Two random 
numbers are generated, u 9 and u p . Then 

<p c = 2%u<p (13) 
20 p c = R^T p (14) 

X c = Xo + pc cos<p c (15a) 
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and 

Y c =Yo + pc sin<p c . (15b) 

Dropping From A User-Specified Xy Position 

5 The condition of dropping a selected particle from a user-specified X-Y 

position is a subset of the previous discussion for the special case of R = 0. The 
positions are stored in the pack surface arrays for that particle until it is placed, at 
which time the drop XY position is replaced by the XY position of the final resting 
location, or the placed particle location, for that particle. 

10 Dropping Above Lowest Sites 

In the preferred implementation, a running list is kept of all placed pack 
spheres whose north poles have not yet been covered by other pack spheres. These 
north poles are ordered according to their altitude (north pole altitude, not sphere 
center altitude). It is generally desirable for the current or selected sphere to find a 

15 final resting place with the lowest altitude. The lowest exposed north pole is a more 
likely candidate for a low final placed particle location. Therefore, if the current 
sphere is to be dropped m times (m > 1), then we choose the lowest m exposed north 
poles and drop randomly within a circle of radius MIN(a c , at) of that north pole's XY 
position, as long as it is within the current sphere's corresponding cylinder wall. 

20 The quantity a c is the radius of the current sphere and ai is the radius of the placed 
sphere above whose exposed north pole the current sphere is being dropped. 
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Dropping randomly within this circle is done in the same manner as is described 
above. North poles outside Wc are ineligible for dropping in this implementation. If 
there are less than m exposed north poles, the other sites are chosen randomly 
within the cylinder. 
5 THE SUBROUTINE PLACE 

The subroutine PLACE in this implementation takes the current particle 
from its initial drop site to the final resting place or placed particle location 
associated with that drop site. The dropping and rolling preferably is done without 
any dynamics. The particles in this implementation do not bounce, nor do they 

10 gain speed or shoot off the edge of another particles with a parabolic trajectory. It 
is as if the particles are being dropped in a highly viscous fluid so that all inertia is 
absent and the particle creeps to its final resting place. 

There are four subroutines in the presently preferred implementation of the 
PLACE subroutine. They are referred to herein as FOBJ1, FOBJ2, FOBJ3 and 

15 STABLE. FOBJ1 finds the first object the descending current sphere encounters. 
FOJB2 finds the second. FOJB3 finds the third. STABLE determines which 
contacted objects the current sphere will stay in contact with and which it will roll 
away from. The term "object" is used instead of "sphere" because one of the objects 
of contact may be the cylinder base or wall. 

20 A number of terms will now be defined for purposes of this discussion before 

describing the operation of the PLACE subroutine. The "roll corridor" is the path of 
descent of a sphere as it descends along its contacts with one or more spheres and/or 
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the cylinder wall. A good share of the effort of subroutine PLACE is to find the 
objects that overlap the roll corridor. An "object" according to this implementation 
can be one of four things that the current sphere can encounter as it descends: (1) 
other pack spheres (placed particles), (2) the cylinder base or wall, (3) a "water 
5 level," or (4) a "catch net." The latter two are described more fully below. "Contact 
stability" refers to whether the current sphere is in compressive contact with an 
object or in tensile contact. If it is in compressive contact with another sphere or 
the cylinder wall, they are pushing against each other due to the weight of the 
current sphere when gravity is acting along the -z axis. The current sphere at this 

10 position stays in contact with an object it compressively touches. Usually if it 
touches three or more objects compressively, it is stable and is placed at that 
position. If there is a tensile contact, the object and current sphere are pulling 
away from each other. Therefore, the current sphere will roll away from the object 
unless it already has three or more compressive contacts. The first pack sphere 

15 touched by the current sphere is denoted as Sphere 1, the second as Sphere 2, and 
so forth. If a set or plurality of spheres are touched simultaneously, the numerical 
order can be arbitrarily selected. 

The flowchart presented in Fig. shows how Subroutine PLACE operates. 

Subroutine FOBJ1 searches for the first object of contact. If the current sphere 

20 simultaneously contacts more than one object, it checks them for stability through 
subroutine STABLE. Subroutine FOBJ2 searches for the second object of contact. 
If the current sphere simultaneously contacts two or more objects, it checks them 



for stability. Likewise, subroutine FO JB3 searches for the third object of contact. 
Again, if the current sphere simultaneously contacts three or more objects, it checks 
them for stability. The current sphere is placed if it is stably touching three or 
more objects or has hit the catch net or a water level. 
5 We will now give an overview of each subroutine and how they interact, 

followed by a detailed discussion of each one. 
Overview of Subroutine STABLE 
Subroutine STABLE determines whether the current sphere is in 
compressive or tensile contact with all objects it is currently touching. It marks 
10 those objects in tensile contact so that the current sphere will roll away from them 
as it continues to descend. 

Overview of Subroutine FOBJ1 
The first subroutine for determining object contact is called FOBJ1 for "Find 
Object 1," object 1 being which is the first object the descending current sphere will 
15 encounter. There are four possible outcomes of FOBJ1; 

• 1. The current sphere encounters the catch net or a water level before 
encountering any other sphere. If it does, then it is placed at the position where 
such contact occurs. 

• 2. The current sphere encounters a placed sphere of the pack. If it does, the 
20 index number of that sphere (hereafter denoted as Sphere 1) is recorded and 

FOBJ2 is called. 
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• 3. The current sphere simultaneously hits two placed spheres of the pack. If it 
compressively touches both of them, it will roll along the roll corridor defined by 
them and FOBJ3 is called to determine whether a third object is hit before it 
loses contact with one or both of the touched placed spheres. 
5 • 4. The current sphere simultaneously hits three or more placed spheres of the 
pack. STABLE checks to see if any three of them provide a stable placement. If 
they do, the current sphere is placed there. If they do not, STABLE marks the 
objects from which the current sphere will roll away. 



10 Overview of Subroutine FOBJ2 

FOBJ2 searches for the second object to be encountered by the current sphere 
as it descends along the path of steepest descent on Sphere 1. There are four 
possible outcomes of FOJB2: 

• 1. The current sphere hits the catch net or a water level before encountering 
15 any other object. If it does, it is placed there. 

• 2. The current sphere loses contact with Sphere 1 by reaching its equator before 
encountering any other object. If it does, it is free falling again and FOBJ1 is 
called, and re-initialized. 

• 3. The current sphere hits a second object, either another placed sphere or its 
20 cylinder wall. STABLE is then called to see if it is in compressive contact with 

both objects. If it is, FOBJ3 is called. If it loses contact with one object, FOBJ2 
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is re-called and re-initialized. If it loses contact with both objects, it is in free 
fall again and F0BJ1 is re-called and re-initialized. 

• 4. The current sphere hits two or more additional objects simultaneously. 
STABLE is called to check for compressive contacts. If there are three or more 

5 compressive contacts, the particle is placed where it is. If there are two, FOBJ3 
is re-called; if there is one, FOBJ2 is re-called; and if there are none, the current 
sphere is in free fall and FOBJ1 is re-called. 
Overview of Subroutine FOBJ3 
FOBJ3 searches for the third object the current sphere encounters as it 
10 descends along the roll corridor formed by the first two objects. There are four 
possible outcomes of FOBJ3. 

• 1. The current sphere hits the catch net or a water level before encountering 
any other object. If it does, it is placed there. 

• 2. The current sphere loses contact with one of its two contacts before 
15 encountering any other object. If it does, FOBJ2 is re-called. 

• 3. The current sphere loses contact with both objects simultaneously before 
encountering any other object. If it does, it is in free fall and FOBJ1 is re-called. 

• 4. The current sphere encounters one or more additional objects. If it does, 
STABLE is called to check whether three or more are in compressive contact. If 

20 they are, the current sphere is placed there. If not, contacts are lost with those 
objects in tensile contact. If two objects remain in tensile contact, FOBJ3 is re- 
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called; if one object, then FOBJ2 is re-called; and if all contact is lost, the current 
sphere is in free fall and FOBJ1 is re-called again. 
Detail of Subroutine STABLE 
The structure and manner of operation of the STABLE subroutine according 
5 to the preferred implementation will now be addressed in greater detail. Consider 
the current sphere with weight W acting in the -z direction. If it is touching one or 
more objects, is it compressively touching them or would it roll away from them. 

Index numbers of placed pack spheres being touched are passed to subroutine 
STABLE in an array called INDEXTOUCH. A separate parameter called 
10 TOUCHWALL states whether the cylinder wall of the current sphere's submode is 
being touched. The manner in which the preferred implementation determines the 
compressive touching criteria when the current sphere touches (a) one sphere, (b) 
two spheres, (c) three spheres, (d) one sphere and the wall, or (e) two spheres and 
the wall will now be addressed. There is generally no need to consider compressive 
15 touching for more than three objects because the case for more than three objects 
can be broken into sets of three. At least three objects should be checked, however, 
because that is the minimum number needed for stable placement in the presently 
preferred implementation. 

Compressive touching of one sphere 
20 Determining compressive touching on one and only one pack sphere can be 

done as follows. If the current sphere is touching its northern hemisphere, it is 
compressive. Because current spheres are dropped from above, they cannot initially 
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touch southern hemispheres and FOBJ1 rejects a sphere as a touching object if the 
current sphere only touches its equator. 

Compressive touching of two spheres 
Consider a current sphere that has descended along a longitude line of 
5 Sphere 1 until it reaches Sphere 2. When it touches Sphere 2, there are only two 
possibilities. Either it will roll along a roll corridor formed by compressively 
touching both Spheres 1 and 2 or it will over-roll Sphere 2 and lose contact with 
Sphere 1. To check for compressive contact with both Spheres 1 and 2, the 
preferred implementation checks whether the current sphere would contact Sphere 

10 1 if it were to roll down the longitude line of Sphere 2 where the current sphere has 
just contacted Sphere 2. If the current sphere would contract Sphere 1, then it 
must be in compressive contact with Sphere 1. It is certainly in compressive contact 
with Sphere 2 because of the contact with it. Therefore, it is in compressive contact 
with both spheres. If, on the other hand, the current sphere would roll away from 

15 Sphere 1 while rolling down the longitude line of Sphere 2, then it will leave Sphere 
1 and is only in compressive contact with Sphere 2. 

Consider the positions of the three spheres: (Xi, Yi 9 Zi) for Sphere 1, (Xs, Y2, 
Zi) for Sphere 2, and (X c , Y c , Z c ) for the current sphere. The radii of these three 
spheres are ai, a>2, and a c , respectively. The Sphere 2 longitudinal line at which the 

20 current sphere is touching Sphere 2 is at azimuth 

(p c2 = tan-i[(7 c - Y 2 )I{X C - Xs)] (16) 
and it touches at latitude 
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9c2 = COS _1 [(Z c - Zs)/(Oc + G2)] 



(17) 



If 8c2 > tt/2, then the current sphere has rolled into the underside or southern 
hemisphere of Sphere 2. In this case, it will be in compressive contact with the two 
5 spheres. If 9 C 2 < n/2, then we must see if it contacts Sphere 1 if it were to roll down 
the longitude line at <p C 2. 

Let the distance between the current sphere and Sphere 1 be denoted D c i 
which is a function of 9 C 2. Its square is given by 

D 2 cX =(X c -X x ) 2 +(7 C -Y x f + (Z C -Z X ) 2 =[(a e +a 2 )sin0 c2 cos^ 2 + X 2 - X x f + 
[(a c + a 2 ) sin# c2 sin^ 2 + Y 2 ■ Y x f + [(a c + a 2 ) cos^ c2 +Z 2 -Z X ] 2 . 

10 (18) 

If Dd or its square increase as G C 2 increases {i.e., as the current sphere rolls down 
Sphere 2 ? s longitude line), then it has lost contact with Sphere 1. That is, if its 
derivative with respect to 9 C 2 is greater than or equal to zero, it is not in 
compressive contact with Sphere 1. If the derivative is negative, it is in 

15 compressive contact with both spheres. The expression below is the derivative to 

within a positive constant of 2(a c +a2), so if it is negative, the current sphere 

compressively touches both Spheres 1 and 2: 

[(a e + a 2 )sm0 c2 cosq> c2 +X 2 - X x ]cos9 c2 cos<p c2 + [(a c + a 2 )sm8 c2 sin^ c2 + 
Y 2 - Y x ]cos0 c2 sin^ c2 -[(a c +a 2 )cos0 c2 + Z 2 -ZJsin^. 

(19) 
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These two-sphere calculations are done in this implementation by Subroutine 

ST2S. 

Compressive Touching Of Three Spheres 
In the case of compressive touching of three spheres, the current sphere at 
5 position (X Cj Y c , Z c ) and with radius a c is in contact with three pack spheres at 

positions (Xi y Yi, Zi), (X2, Y2, Z2), and (X3, Y3, Z3) and with radii ai f 0,2, and a?. Let 
the current sphere have weight Wdue to gravity or similar force. Assume it is 
attached to the three spheres at the point of contact with them. The preferred 
implementation solves a force equation to determine whether the attachments are 
10 compressive or tensile. If an attachment is tensile, the current sphere will roll 
away from it. 

Let the unit vectors pointing from the current sphere toward Sphere i where i 
is 1, 2, or 3, be denoted by eu Let the magnitude of the force applied on the current 
sphere by Sphere i be denoted by fi. If fi is positive, the attachment is tensile; if 
15 negative, it is compressive. The three equations needed to solve for the three fi are 
given by the vector equation 

fiei + fseM + fses = (0,0 -W) . (20) 
Because in the implementation we are only interested in whether the fs are 
compressive (fi < 0) or tensile (fi > 0) and are not interested in their magnitudes, the 
20 preferred implementation simply sets Wto unity and obtains for the / solutions 

fi - (e3xe2y - e2xe3y)/ det E 
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f2 = (eix e3y - e3x ei y )/det E (21) 
fs = (€2xei y - eixe2y)f det E 

where det E is the determinant of the 3x3 matrix whose columns are the unit 
vectors ei, e2, and es: 
5 det E = eix(e2y esz - e3y e2z) + e2x (esy eiz - ei y e3z) + e3% (eiy e2z - e2y eiz) * (22) 

If the determinant is zero, the set of equations has no solution because the 
three placed pack spheres and the current sphere all lie in the same plane. In the 
unlikely event that this should happen in practice, the preferred implementation 
first sees if another set of three spheres gives stability if the current sphere is 
10 touching more than three placed pack spheres. If there is no stable triplet, then it 
enacts a subroutine to perturb the current sphere's position slightly in a downward 
direction such that it will not overlap any of the three spheres. If the perturbation 
causes contact to be lost with one or two of the touching spheres, it can proceed with 
FOJB2 or FOJB3. If no contact is lost, then at least the four spheres are no longer 
15 coplanar and it can redo the processing without a zero determinant. 

These three-sphere calculations are done in the preferred implementation by 
Subroutine ST3S. 

Compressive Touching Of One Sphere And The Wall 
In the case where there is compressive touching of one sphere and the wall, 
20 the wall cannot be the first object contacted because in this implementation, the 
current sphere drops parallel to it when free falling. Therefore, the only way the 
current sphere can be in contact with both its cylinder wall and one placed pack 
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sphere is if it hits Sphere 1 and subsequently rolls into the wall. The current 
sphere is always in compressive contact when rolling on one sphere. Furthermore, 
because the current sphere must be in contact with Sphere l's northern hemisphere 
when it contacts the wall ? it must be in compressive contact with the wall because it 
5 just rolled into it. Therefore, any time the current sphere is in contact with one 
sphere and the wall, it is always in compressive contact with both and will descend 
along a roll corridor defined by both. In the unlikely event that the current sphere 
hits the cylinder wall and nothing else at the same time it reaches the equator of 
Sphere 1, it free falls from that point and FOBJ1 is re-called. 

10 Compressive Touching Of Two Spheres And The Wall 

The case of compressive touching of two spheres and the wall is a simple 
offshoot of the three-sphere case. The cylinder wall will replace one of the spheres 
as a contact. However, in the absence of friction, the wall can only apply a force to 
the current sphere directed radially inward. 

15 Let ri and r2 denote the positions of the two pack spheres and r c denote the 

position of the current sphere. The relative positions of the two pack spheres to the 
current sphere is n c = n - r c , i = 1 or 2. These relative vectors point from the 
current sphere center to the i-th pack sphere center. Let etc be the unit vector 
pointing along n e . The force the current sphere applies to each pack sphere and the 

20 wall must sum to its weight W. That gives the equation 

fieic + fzesc + fwewc = (0,0 -W) . (23) 
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The unit vector ewe, which points radially outward from the cylinder center, is given 
by 

ewc = ( Xc i + rcMxc 2 + yc 2 ) 1 ' 2 . (24) 
As before in the three sphere case, the preferred implementation solves for the fi, i 
5 = 1, 2, or w, but actually only their sign is needed: positive, negative or zero. Hence 
Wis set to unity. Let eix be the x component of the vector etc, i - 1, 2, or w, and 
similarly for the y and z components. Then the solutions for the fi are given by 
fi — (e2y e W x - e2% e^)/det e 

f2 = (ei x e W y - ei y e W x)f det e (25) 
10 f w = (eix e2y - ei y e2x)/det e 
where 

det e = -eix eu;y + ei^ e^x + eiz(e2x e W y - e2y e W x) . (26) 
If an fi, i=l,2, or w, is negative, the force between the current sphere and that 
object is compressive and they stay in contact. If it is positive, the force is tensile 
15 and the current sphere will roll away from that object. In the unlikely event that it 
is zero, the implementation considers it tensile and has the current sphere roll 
away from the object, thereby trying to find a point of stability at lower altitude. If 
all three fi are negative, the sphere is in stable contact and is placed at its current 
position. 

20 These stability calculations for two spheres and a wall are performed in the 

preferred implementation by Subroutine ST2S1W. 
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Detail of Subroutine F0BJ1 
Returning to the subroutine FOBJ1, it finds the first object the descending 
current sphere encounters as it is "dropped" from the position selected by 
subroutine DROPSITE. There are three possible objects the current or selected 
5 sphere can encounter: (1) an already-placed pack sphere, (2) a water level, or (3) 
the catch net. A fourth object that may be encountered later is the cylinder wall, 
but because the current sphere is assumed to drop parallel to the cylinder wall, it 
cannot contact the wall until it rolls off some placed pack sphere. 
The Catch Net 

10 In accordance with the method aspect of the invention, the selected particle 

placement includes defining a catch net representative of buoyancy of a portion 
of the placed particles and positioning the catch net within the space based upon 
the placement of the portion of the placed particles. The selected particle is placed 
in non-overlapping relation with respect to the catch net. Preferred versions of the 

15 method include defining a pack surface for the placed particle, and positioning the 
catch net relative to the pack surface. The pack surface may comprise a top layer of 
the placed particles within the space, or preferably, within each of the subspaces. 
The particle surface optionally may correspond to an average of the south poles or 
south pole positions of the top layer of placed particles. The catch net may be 

20 positioned at a distance from the pack surface based upon a selected particle radius, 
such as an average particle radius of the top layer or another subset of the placed 
particles at the pack surface. 
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In the way of background, when simulating a particle pack, particularly one 
with a wide particle size distribution, if the volume fraction of large particles is 
much larger than that of small particles, the small particles can rattle down 
through the pack by falling through the interstitial spaces of the large particles. In 
5 many particulate materials, there are buoyancy or suspension considerations that 
would keep the small particles more or less well distributed throughout the pack. 
Examples of such buoyancy or suspension factors may include the nature or 
properties of the binder or other inter-particle medium or media, the interaction of 
the particles with the inter-particle medium or media, interactions among the 
10 particles themselves (such as electrical or magnetic interactions), interactions with 
externals forces or fields that tend to disperse or suspend the particles, and similar 
phenomena. 

In accordance with this aspect of the present invention, this buoyancy or 
suspension phenomena is simulated using a "catch net." The catch net comprises 
15 an assumed surface or construct to catch the current or selected particle as it is 
placed and to prevent that particle from passing a defined location or level which 
theoretically represents space occupied by smaller ones of the previously placed 
particles. 

In the presently preferred implementation, a catch net is incorporated into 
20 the space, below which the south pole of the falling particles cannot descend. When 
a particle's south pole hits the catch net belonging to its cylinder, it is placed in the 
pack at that point, even if it is not touching any other particle. Preferably there is a 
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separate catch net for each particle submode, and thus for each subspace or 
cylinder. Each catch net preferably covers the entire cross section of that submode's 
cylinder. If the current sphere's south pole hits that catch net, it is placed there. It 
does not see or respond to any other submode's catch net. The altitude of the catch 
5 net preferably is determined dynamically, that is, it can depend on which placed 
particles are being contacted. 

The use of a catch net is not necessarily limited to a specific type or size of 
particle, or a particular size or geometry of the space or subspaces. In accordance 

10 with the presently preferred implementation, however, as noted, the particles 
comprise spheres and the characteristic dimension of each of the particles 
comprises a radius. The central string comprises a line, preferably the z axis, 
within the space. The space comprises a cylinder, and the subspaces comprise 
concentric cylinders disposed about the central string, each of the cylinders 

15 defining one of the subspaces and having a radius that corresponds to the radius of 
the particles in the category corresponding to that cylinder. 

In accordance with the presently preferred implementation, the method 
further includes defining a pack surface for the placed particles, and the catch net 
positioning comprises positioning the catch net relative to the pack surface. Also in 

20 accordance with this implementation, the catch net positioning comprises 
positioning the catch net a distance from the pack surface based upon a 
selected particle radius. 



The placed particles may be assumed to comprise a top layer, and the pack 
surface optionally but preferably comprises an average of the particle locations of 
the placed particles in the top layer. In this implementation, the particle surface 
corresponds to the south poles of the top layer placed particles, and more 
5 specifically corresponds to an average of the south pole positions of the top 
layer placed particles. 

In the presently preferred implementation, the catch net extends across the 
cross sectional area of the space, across the cross sectional areas of each of the 
subspaces. 

10 In this preferred implementation, the catch net comprises N subnets, one of 

the subnets corresponding to each of the subspaces. Each of the concentric 
cylinders 1, 2, 3, .. N has a circular cross section perpendicular to the central string. 
In the preferred implementation, the catch net comprises a plurality N of subnets, 
or individual catch nets, one for each of the cylinders. Each of the subnets 

15 preferably extends over the cross sectional area of the corresponding subspace or 

cylinder. The position or level of the catch net or subnet for each of these subspaces 
can and typically will differ from one another, although in some instances they may 
be at the same or nearly the same level. 

Preferably, the positioning of the catch net positioning comprises setting the 

20 catch net or each of the subnets at a selected distance from the top surface. This 
selected distance may be a user-defined quantity, or it may be dictated by the 
material being simulated or other factors. The catch net preferably is positioned at 
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an altitude along the z-axis that is a function of the particle radius of particles 
corresponding to the mode or submode of the cylinder in which the catch net or 
subnet resides. The catch net or subnet preferably but optionally is some multiple 
of the particle radius, preferably one or two ratio, below the top layer or pack 
5 surface. Adjustments or offsets, however, may be used. 



When the pack is first being formed so that not enough particles have yet 
been placed to adequately define a pack surface, the catch net may be defined as the 

10 bottom of the cylinder into which the particles are being dropped. In some 

instances, however, particularly where the specific application involves simulating 
a very large three-dimensional pack, it is often advantageous or desirable to avoid 
the edge effects of the flat bottom of the cylinder. Accordingly, in some instances 
the positioning of the catch net may comprises spacing the catch net from the base 

15 surface of the space. The spacing from the base surface, for example, may simulate 
the positioning of the catch net for a top layer of the placed particles, for 
example, the positioning of the catch net off of the base of the space may be done in 
a way that the bottom layer of placed particles simulate the shape, configuration 
and appearance of the top layer of particles after the particle pack has been 

20 constructed to include multiple layers. 

This may be implemented by assigning as initial catch net position for a 
kth one of the subspaces Z in it(k), assigning as the characteristic dimension of 



the particles of a kth one of the particle categories ak, assigning as the 
characteristic dimension of a small one of the particles amin of the particles, 
and assigning as the characteristic dimension of a large one of the particles 
amax. Below a threshold level the catch net positioning further comprises 
5 positioning the catch net for a kth one of the particle categories at a catch net 
position Znet(k) within a kth one of the subspaces determined by 

ZnetQt) — Zinit + H V CLk (Xminl Clmax 

where r represents a weighting coefficient and H represents a switching 
coefficient. The weighting coefficient preferably but optionally is assigned a 
10 random number. Below the threshold value the switching coefficient preferably is 
assigned a value of one, and above the threshold value the switching coefficient 
preferably is assigned a value of zero. 

In the presently preferred implementation, the switching coefficient 
mechanism is implemented by a switch called ROUGHNET which, when true, 
15 invokes the random number generator to slightly change the altitude of the catch 
net so that the first set or layer of placed particles do not all have the same altitude. 
If ROUGHNET is false, the south poles of the first set of placed particles do have 
the same altitude because they hit the same catch net. 

Let Zinit be the initial altitude above which the pack will be built. It is 
20 specified by the user in the preferred implementation but its default preferably is 
zero. The initial catch net altitude for the &-th mode is 

ZnetQt) ~ Zinit + H r dk CLminlamax (30) 
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where r is a random number in [0,1], a m in is the smallest particle radius of all 
submodes belonging to the pack, Umax IS the largest, and His unity if the 
ROUGHNET switch is true and zero if it is false. This expression for the catch net 
altitude holds until enough particles have been placed to establish a pack surface, 
5 i.e., until m (which is greater than or equal to m m in) is reached. No part of any 
particle is permitted to descend below Zinit in this implementation. 

The catch net preferably is positioned relative to the pack surface. The pack 
surface can be determined in a number of ways. It may be determined, for 
10 example, based on the large particles in the top layer, or another subset of 
the placed particles, e.g., such as the last particles placed. 

The particle surface may be ascertained in the following manner. For each 
particle category k of the placed particles in the top layer, the program defines a 
particle radius ai for the placed particles i of that category k. For the cylinder k 
15 corresponding to the particle category k, the programs assigns a cylinder radius 
Wk. The program also assigns a top layer particle number m(k) and determines 
values for m(k) by evaluating 



where N is the number of particle categories, e.g., the number of submodes. The 
20 program then determines the pack surface location using 




Submode(i)<k 



Submode(i)<k 
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1 m 

where S represents the pack surface and Zi represents the position of a center of a 
center one of the placed spheres. 

In accordance with the preferred implementation, the surface altitude S of 
5 the pack is defined as the average south pole altitude of the last m placed particles 
where m is a positive integer that will be discussed below. This surface altitude is 
given by 

tYl i=l 

where Zi is the altitude of the center of Sphere i which is one of the last m placed 
10 particles and at is its radius. Instead of using a single ith particle radius in 

Equivalence. 27, one may generalize the expression by permitting (Zi - ai) to be 
written (Zi - f(ai)) where f(a0 is a function of particle radius ai. The preferred value 
for f(ai) = ai, but f(a0 may assume other forms, for example, such as f(ai) = nai where 
n is a positive real or integer number. 
15 Although the "surface" of a random particle pack may be a defined in a 

number of ways, the surface preferably is defined only after the integer m is large 
enough that the falling spheres have covered the cylinder area. The A;-th cylinder 
with radius Wk has area nWk 2 . Let m be defined as the number of most recently 
placed spheres whose submodes have cylinder radii less than or equal to Wk and 
20 whose summed cross-sectional area just exceeds or equals the &-th cylinder area. 
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For the k-ih submode, m is taken from the following condition. Let m(k) be the 
number m for the k-ih submode. It is found from 



where iVis the total number of submodes. The submode k for which this condition 
5 is satisfied first, that is, has the lowest m(k) becomes the m for the pack: 



In the preferred implementation there is a minimum number mmin for those 
unused packs that defeat the reasoning leading to the value for m given in Eqs (28) 
and (29). In the cases the default value of mmin to 12 is assigned in the preferred 
10 implementation. 

Once a surface is established, its altitude S is determined by using Eq. (27). 
As the current sphere descends, if it hits a pack sphere much larger than itself, and 
if the TUNNEL switch is true, it will tunnel through the large sphere to check for 
possible cavities under the large particle. Therefore the catch net should be placed 
15 quite low so the current sphere can investigate possible cavities. If TUNNEL is 
false, the catch net does not need to be that low. 

If the current sphere hits a sphere somewhat bigger than itself, then it is 
unlikely to tunnel because the surface sphere is not large enough to accommodate it 
in a cavity. Then the catch net need not be as deep. 



Submode(i)<k Submoie(i)<k 




(28) 



m = MIN[m(fe)] . 



(29) 
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Finally, if the current sphere hits a sphere smaller than itself, then it cannot 
penetrate too deeply into the pack so the catch net can be rather shallow. 

Letting the subscript i denote the first pack sphere contacted, the catch net 
criteria are quantitatively summarized as follows: 
5 If aila c < 1, then the catch net is shallow: 

Znet(k) = S-ai. (31a) 
If 1 < ailac < a x , where a x preferably is V6 + 2, then the catch net is moderate: 
Znet(k) = S-2a c . (31b) 
If ailac > a x , where ax preferably is V6 + 2, then the catch net is deep: 
10 Znet(k) = S-2a c -ai. (31c) 

The value of a x represents the sizes required for a current sphere to fit into a cavity 
formed by larger spheres, which preferably assumes a valueof(V6 + 2). If the 
centers of four identical larger spheres form a regular tetrahedron, and a smaller 
sphere just fits into the interstitial space of the tetrahedron, then the ratio of the 
15 smaller sphere to larger sphere diameter is V6 + 2 in preferred implementation. 

Further in accordance with the preferred implementation, the code monitors 
the selected or current particle placement to identify the occurrence of the following 
sequence of events: 

1. The catch net is set at an altitude when the current sphere hits the i-th pack 
20 sphere. 

2. Then the current sphere rolls off the i-th sphere and loses contact with it. 
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3. When it hits another pack sphere, the catch net is set above the current 
altitude of the current sphere's south pole. 

If this series of events is detected, the current sphere is raised (its position is 
adjusted upward) until its south pole is at the new catch net altitude, and it is then 
5 placed there. 

The Water Level 

In accordance with another aspect of the invention, the placement of the 
particles comprises defining a water level representative of a level of a portion 

10 of the placed particles that are smaller than the selected particle and 

represent a surface of the smaller placed particles, and positioning the water level 
within the space based upon the smaller particle surface. The current or selected 
particle then is placed in non-overlapping relation with respect to the water level. 
The water level serves a purpose different than the catch net. Recall that 

15 the particle pack according to the preferred implementation is constructed using a 
series of concentric cylinders, one cylinder defined for each particle category or 
submode. Recall further that the centers of particles in this implementation cannot 
pass outside their own corresponding and constraining cylinder into which they are 
dropped. Large cylinders are used for large particles and small cylinders for small 

20 particles, unless the user for some reason designates otherwise. The cylinder 

construct facilitates or permits use of the reduced dimension approach to particle 
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pack construction. Its purpose is to prevent the program code from having to fill 
large cylinder volumes with potentially very large numbers of small particles. 

This approach, however, creates the problem that the large cylinders will 
contain empty space into which large spheres can roll during simulated particle 
5 placement where, in a corresponding physical particle pack such as the one being 
simulated, this would not occur. In the actual pack, the smaller particles would be 
included even in the regions comprising the larger cylinders, and these smaller 
particles would be encountered to prevent such rolling. 

The water level simulates this presence in the larger cylinder areas of 
10 smaller particles that, according to the presently preferred implementation, have 
not been placed in the larger cylinders. 

There are a number of methods or approaches for selecting or setting the 
water level. According to one such approach, the water level positioning comprises 
determining an average location of the particle locations along the central string 
15 of the particles of the portion of placed particles that are smaller than the 

current or selected particle, and positioning the water level at the average location. 
The water level may be set at the average altitude of the small particle 
"surface" as if the outside cylinders were filled with small particles. Then when a 
large particle or sphere descends onto the water level, the particle, and preferably 
20 its south pole, is stopped there as it would be if that volume had been filled with the 
smaller particles. 
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In the presently preferred implementation, wherein the plurality of 
subspaces (here concentric cylinders) are used, the water level comprises a plurality 
of subspace water levels, preferably one subspace water level for each subspace or 
cylinder. These water level positions may and typically will differ from one another 
5 as the particle pack is constructed. The water level positioning optionally but 

preferably comprises assigning a subspace water level position to a portion of 
the subspace water levels. The portion of the subspace water levels preferably 
correspond to one water level for each cylinder except the first cylinder. In the 
preferred implementation, all of the cylinders except the first or central one, 
10 containing the smallest particles, has a subspace water level and corresponding 
subspace water level position. 

In the presently preferred implementations, the water level positioning 
comprises using an offset to position the water level, or an offset relative to the 
average location of the particle locations. Preferably, the water level positioning 
15 comprises using an offset for each of the subspace water level positions. 

The water level positioning similarly may comprise using an offset for each of the 
subspace water level positions. 

In the presently preferred implementation, the water level altitude for the i- 
th cylinder (i > 1) is defined as the average exposed north pole altitude of all 
20 particles with radius less than au The user may raise or lower this level by adding 
a user-defined offset, but the default for this offset in this implementation is zero. 
This water level exists between the £-th cylinder wall of radius Wi and the next 
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smaller cylinder wall of radius WLi. Thus the spaces between different cylinder 
walls have different water levels. The cylinder corresponding to the smallest 
particle's submode has no water level, only a catch net. 

If a larger sphere is descending onto the pack near the center cylinder, it may 
hit the center cylinder's pack of small particles and tend to roll off its edge and fall 
to the floor of the pack or fall until it hits another large sphere. Thus the pack 
would segregate with larger particles outside smaller particles. With the water 
level in place, however, if the larger sphere rolls off the pack in the central cylinder, 
it simply "floats" on the water level which is at the same level as the "surface" of the 
central cylinder pack. Thus, any time a descending sphere hits a water level, it is 
placed where it is, just as was done with the catch net. Ordinarily, water levels will 
be higher than catch nets, although this is not always true. 

Search for the First Object Encountered 

A free-falling current sphere can encounter the following situations. It might 
hit a single pack sphere. If so, FOJB2 is called. It can hit two pack spheres 
simultaneously. If so, FOBJ3 is called. It can hit three or more pack spheres 
simultaneously. If so STABLE is called to see if it is in stable contact with any 
triplet of these spheres. Hitting two or more pack spheres simultaneously is highly 
improbable in a double precision random code but not so improbable if the pack is 
semi-ordered, and hence these scenarios are included in the preferred 
implementation. 
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In addition, the south pole of the current sphere may hit the catch net or a 
water level before it hits any other sphere. In either case, it is placed where it is. 

A pack sphere, e.g., they-th sphere, is a candidate for being the first object 
encountered by the descending current sphere if their two radii overlap, i.e., if 
(Xc - Xj)2 + (Y c - < (etc + a;)* (32) 
If there were no other obstacles in the way, the current sphere would strike Sphere j 
at an altitude of 

Zc(j) = Zj + [(oc + a;)* - (Xc - Xj) 2 - (Y c - W 2 . (33) 

Once the pack has grown to many spheres, efficiency can be greatly enhanced 
to the extent that as many spheres as possible be excluded from the detailed 
considerations of the above two equations, while not unduly adversely affecting the 
accuracy of the simulation in representing the physical particle pack being 
simulated. For example, if | x c - x; I or \y c -yj\ alone were greater than (a c + oj), 
then there is no need to check the square of their sums as above. 

The processing demands are reduced in the preferred implementation 
because all spheres whose north poles have submerged beneath the catch net are off 
limits and need not be included. Simply flagging them by putting a minus sign in 
front of their mode or submode number indicates that no calculations are to be done 
to see if they are in range of the descending current sphere. Furthermore, in this 
implementation the program keeps track of the lowest index number for any sphere 
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whose north pole protrudes above that mode's catch net so that the eventually large 
list of particles beneath the catch net are not searched at all. 

Submergence of a placed pack sphere beneath the water level, however, is not 
necessarily grounds for dismissing it from consideration. The current sphere may 
5 roll into a cylinder for which the water level is different and the pack sphere may 
again become a candidate for contact. 

If the current sphere is dropped far enough outside a set of cylinders that 
there is no chance for a first contact with the spheres belonging to those cylinders, 
then in the preferred implementation they are dismissed from further 

10 consideration. If the submode number of a sphere belongs to an out-of-range 

cylinder, no further processing is required in the preferred implementation. Let Wk 
be the cylinder wall radius for the k-th particle submode. If 
(Wk + ak) 2 <Xc 2 +Yc 2 (34) 
then no fe-th mode sphere is a candidate for first contact with the descending 

15 current sphere in this implementation. 

In summary then, FOBJ1 searches the list of all placed pack spheres above 
the catch net for the submode to which the current sphere belongs, excluding all 
placed spheres that may be in interior cylinders that the current sphere cannot 
contact or overlap. The possibility of overlap in the x directions is checked, followed 

20 by the possibility of overlap in the y directions, followed by a direct calculation of 
overlap. The altitude z c (j) at which the current sphere would contact Sphere j is 
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stored and that sphere is rejected if other placed pack spheres have a higher contact 
altitude. Again, if no pack spheres make contact before the current sphere reaches 
the catch net or water level, it stops and is placed there. 
Detail of Subroutine FOBJ2 
5 Subroutine FOBJ2 finds the second object to be contacted by the descending 

current sphere. FOBJ2 would only have been called if the current sphere had first 
contacted a single pack sphere that has been denoted Sphere 1. If it had contacted 
the catch net or water level, the current sphere would have been placed there and 
FOBJ2 would never have been called. If it had contacted more than one sphere 

10 simultaneously, FOBJ3 and/or STABLE would have been called. FOBJ2 is called 
only if a single pack sphere's northern hemisphere has been touched. 

When the current sphere touches the northern hemisphere of Sphere 1, it 
begins to roll down the path of steepest descent which is Sphere l's longitudinal 
line passing through the point of contact. The volume swept out as the 

15 current sphere descends along this path toward the equator is called the roll 

corridor. The end of the roll corridor is either where the current sphere's center 
reaches the equatorial plane of the Sphere 1, reaches its cylinder wall, or reaches 
the catch net or water level, whichever is higher. The preferred implementation 
identifies all other pack spheres that protrude into the roll corridor. 

20 Processing can be somewhat simplified if the cylinder wall is assumed to be 

the boundary that the center of the current sphere cannot pass, rather than the 
boundary through which its leading edge equator cannot pass. 



If the center of the current sphere intersects the cylinder wall during descent, 
then that point is the end of the roll corridor because it will at least have found the 
wall as its second object if no other pack sphere has been contacted first. 

The beginning of the roll corridor is then the place of initial contact between 
5 the current sphere and Sphere 1. The current sphere's polar angle at that initial 
contact is given by 

9 ci = cos _1 [(^c - Zi)/(ac + ai)] . (35) 

The end of the roll corridor 0 c f will either be at the equatorial plane of Sphere 
1 or where the center of the current sphere intersects the catch net, the water level, 
10 or the cylinder wall. If it is at the equatorial plane, then 

Gcf=7u/2. (36) 
If the altitude of the catch net Z m t > Zi-a c , then 

G c f = COS^KZnet + Gc - Zi)/(a c + ai)] . (37) 
If the altitude of the water level Z w > Zi - a c > then 
15 Gcf = cos" 1 [(Zw + a c - Zi)/(ac + ai)] . (38) 

The last possibility is that the roll corridor ends with an intersection of the 
cylinder wall. Let <p c be the azimuthal position at which the current sphere has 
touched Sphere 1. It is given by 

cp c = tan-i[( Y c - Yi)/(Xc - Xi)] . (39) 
20 The radial position of the current sphere as it descends Sphere 1 as a function of G c 
is given by 
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{[Xi + (a c + ai) sin9 c cos<p c ] 2 + [Yi + (a c + ai) sin9 c sincpc] 2 } 1 ' 2 . (40) 
If this radial position is equal to the cylinder radius Wfo? some value of 9 C between 
0 and 71/2, then the roll corridor intersects the cylinder wall. Setting this expression 
equal to W> squaring both sides, and solving the quadratic equation for sin0 c gives 

. a - ( X x cos <p c + Y x sin <p c ) ± 4{X X cos <p e + Y x sin <p c f +W 1 - X\ - Y* 
5 sin 9 cf s ( 41 ) 

Because 0 C f is a polar angle, it lies in the range [0 ? 7i]. Therefore, sin0 C f can never be 
negative. If both right hand side solutions of Eq. (41) are greater than unity, there 
is no solution and the current sphere does not intersect the cylinder wall. If both 
are less than or equal to unity, then we choose the larger one. 

10 The presently preferred implementation conducts appropriate processing to 

evaluate each of these possible circumstances. Before evaluating Eq. (41), the 
program code according to this implementation adds the global radial position of 
Sphere 1 to the sum of the radii of Sphere 1 and the current sphere. If this sum is 
less than or equal to W 9 then the sphere cannot hit the cylinder wall in any case so 

15 the more involved evaluation above is unnecessary and preferably is not performed. 

The end of the roll corridor is the minimum 0 c f of Eqs. (36, 37, 38, and 41). If 
the current sphere encounters the catch net or water level, it is placed there. 

The preferred implementation also checks for all pack spheres whose volumes 
intrude upon the roll corridor swept out by the current sphere which covers the arc 

20 from 9ci to 0 c f and azimuth <p c on Sphere 1. If there is a 0 C with smaller value than 
the end of the roll corridor 0 c f at which the current sphere encounters one or more 
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pack spheres, then the sphere or spheres touched are noted. Program control is 
then passed to subroutines that deal with simultaneous contact of two or more 
spheres, as described herein. 

If no pack spheres overlap the roll corridor, then the current sphere is placed 
5 if it hits the catch net or water level. If it hits neither of these, then it rolls to the 
equator and drops from there. At that point it is a free falling sphere again so 
FOBJ1 is re-called. 

We now turn our attention to finding which placed pack spheres overlap the 
roll corridor and at what point the current sphere would strike them as it descends 
10 from 0 C i. The term "global coordinate system" is used herein when discussing a 
single coordinate system in which the entire pack is defined. Capital letters are 
used when referring to a global coordinate system, e.g., J(. = {X i9 Y i9 Z t ) for the 
Cartesian position of the i-th sphere. 

A local coordinate system, on the other hand, is one whose origin is at the 
15 center of one of the spheres and whose axes are parallel to the global axes. We use 
small letters to define it, e.g., )r. = (x^y^z^ for a position relative to the center of the 

i-th sphere. 

Consider the j-th pack sphere. Its center is at global position Rj = (Xj, Y h Zj) 
and its radius is aj. In the global coordinate system, the position of the current 
20 sphere as it descends along Sphere i is 

J? c = J(. + (a c + <Z;)(sin 0 c cos q> c \ + sin 9 C sin <p c j + cos 0 c k) . (42) 
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The current sphere rolls into Sphere j if there is a solution 9 C that lies in the 
interval 0 C i < 0 C < 0 c f to the equation 

\i-^\ = a c +aj. (43) 

If this equation can be satisfied at all, it will have two solutions: one where the 
5 current sphere touches Sphere j from above and another where the current sphere 
touches Sphere j from below. (If it just nicks Sphere j, these two solutions can be at 
the same 0 C> but if Sphere j is just nicked, the preferred implementation considers 
that a miss.) We are only interested in the top or smaller 0 C solution, and then only 
if the top 9c lies in the roll corridor range. 

10 Explicitly, the equation to be solved is 

[Xi -Xj + (a c + ai) sin0 c coscp c ] 2 +[Yi-Yj + (a c + ai) sin9 c sincpj 2 + 
[Zi ~Zj + (a c + ai) cosGc] 2 = (a c + aj) 2 (44) 
which, defining Xij = Xi~ Xj y Yij=Yi- Yj, and Zij = Zi- Z h can be rewritten 
(Xij coscpc + Yij sin<p c ) sin0 c + Zij cos0 c = [(aj - at)(aj + at + 2a c ) - Rij 2 ]/[2(a c + at)] 

15 (45) 
where 

Ro* = Xjj* + Yij* + Zij*. (46) 
This equation has the form 

A sin0 c + B cosOc = C (47) 
20 where 

A = Xij coscpc + Yij sincpc (48) 

76 



B = Zij (49) 
C = [(aj - ai)(aj +ai+ 2a c ) - i?y 2 ]/[2(a c + ai)] . (50) 
Dividing through by V(A 2 + B 2 ), we obtain 

cosa sinGc + sina cos9 c = CN( A 2 + B 2 ) (5 1) 

5 where 

cosa = A/V( A 2 + B 2 ) (52) 
and 

sina = BN( A 2 + B 2 ) . (53) 
A common trigonometric identity turns Eq. (51) into 
10 sin(9c + a) = CN( A 2 + B 2 ) (54) 
which has solution 

9c = sin"i [C/V( A 2 + B 2 )] - a (55) 
where 

a = tan _1 [A/B] . (56) 
15 If | CH( A 2 +B 2 ) \ > 1, then there is no solution and the j-th sphere does not overlap 
the roll corridor. 

The FORTRAN function ATAN2(A,£) or its equivalent in other programming 
languages uses information on the signs of A and B to determine in which quadrant 
a lies so it is uniquely defined. The arcsine, however, has two solutions. One 
20 solution is the principal value which lies in [-7r/2,7i/2] and is found by the standard 
FORTRAN function ASIN or its equivalent. The other solution is n minus this 
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principal value. These two solutions comprise the two solutions for the current 
sphere contacting Sphere j on either side along the longitude line defined by 
azimuthal angle cp c . Denoting the principal value by capitalizing the arcsine, the 
two solutions are 

5 e c i = SixTi[C/V( A 2 + B 2 )] - a (57a) 
and 

0c2 = n - SiiTi [CH( A 2 + B 2 )] - a . (57b) 

If either of these 0 lie outside [0,27i), then 2n must be added or subtracted enough 

times to bring the solutions back into [0,27i). If neither of these solutions lie in the 
10 open interval (0 C i, 0 C f), then the current sphere does not contact Sphere j. If one 

solution does lie in (0 C i, 9 C f), that is the position at which the second contact is made. 

If both do, the lesser of the two solutions is the position of second contact. 

Finally, whichever pack sphere has the smallest 0 C is Sphere 2. That point of 

contact is denoted 0 C j. If two or more pack spheres are contacted simultaneously, 
15 then all their index numbers are returned by FOBJ2. If the cylinder wall and one 

or more pack spheres are contacted simultaneously, their index numbers are 

returned and the flag indicating a wall contact is activated. 

As the preferred implementation searches through the pack spheres to find 

which is touched first, there are several efficiency criteria can be used to eliminate 
20 most spheres from the detailed consideration of the above equations. 
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First, if the roll corridor is far enough away from the global Z axis, then the 
current sphere may not be able to approach some inner cylinders closely enough to 
touch their spheres. Therefore, those pack spheres whose mode number indicates 
they belong to the inner cylinders need not be considered further. The preferred 
5 implementation can determine which cylinder's spheres can never intersect the 
current sphere's roll corridor? 

Because the current sphere rolls down a longitude line on Sphere i y when 
viewed from above, its projection is a straight line segment on the XY plane with 
beginning point (Xd, Yd) and ending point (Xcf, Y C f) where 

10 Xd -Xi + (a c + ai) sinGci coscp c (58a) 
Yd - Yi + (a c + at) sinOci sincpc (58b) 
Xcf = Xi + (a c + ai) sinGcf coscpc (58c) 
Y c f = Yi + (a c + ai) sinO c f sincpc . (58d) 
To find whether this line segment comes close enough to the Z-th cylinder to touch 

15 any of its spheres, the preferred implementation turns the current sphere's position 
into a pair of parametric equations: 



where s is a dimensionless parameter that is zero at the beginning of the roll 
20 corridor and unity at the end. In order for the current sphere to have any chance of 
touching an Z-th mode sphere, it must come closer than the sum of its radius plus an 



Xc — Xd + (Xcf— Xd)s 



(59a) 



Yc = Yd + (Ycf- Yd)s 



(59b) 
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Z-th sphere's radius from the Z-th cylinder wall which has radius Wu In other words, 
some portion of the roll corridor must come within Wi+ai+ a c of the global Z axis. 
This condition is expressed by 

Xc 2 + Yc 2 = (Wi + at + ac) 2 (60) 
5 for a value s somewhere in the open interval (0,1). If Eq. (60) is never satisfied, 
then the current sphere cannot touch any of the Z-th mode spheres. Explicitly 
written, Eq. (60) is a quadratic equation in s: 

As 2 + 2Bs + C = 0 (61) 
where 

10 A = (Xcf- Xci) 2 + (Ycf- Yd) 2 (62) 
B = Xci(Xcf- Xci) + Yci(Yc f - Yd) (63) 
C = Xci 2 + Y*?- (Wi +ai+ a c ) 2 . (64) 
The solutions are 

~B±^B 2 -AC 

s = (65) 

15 If B 2 ~AC is negative, there is no solution. If it is zero, there is only one real solution 
which means the roll corridor just grazes the Z-th cylinder, in which case it is 
considered a miss. If it is positive, there are two real solutions. The preferred 
implementation solves for both of them, and checks to determine whether either one 
lies in the open interval (0,1). If one or both lie in this interval, the Z-th mode 

20 spheres are candidates for touching the descending current sphere. 
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To find which cylinders are inaccessible to the current sphere, the preferred 
implementation starts with the innermost cylinder (/ = 1) and continues outwardly 
until it reaches a cylinder that is big enough that it is accessible. All larger 
cylinders also will be accessible. 
5 There also is a simple altitude criteria to delimit detailed consideration of 

pack spheres as candidates for touching the current sphere. The global altitudes of 
the beginning and end of the roll corridor are 

Zd = Zi + (ac + at) cosGd (66a) 
Zcf = Zi + (a c + at) cos9 c f , (66b) 
10 respectively. Therefore, if a j-th pack sphere lies below the end of the roll corridor 
by the amount 

Zj < Zcf- a c ~ aj (67a) 
or lies above the top of the roll corridor by the amount 

Zj^Zci+ac+aj, (67b) 
15 then it is not a candidate for touching the current sphere. 

Likewise, spheres lying beyond the roll corridor in the x and y directions need 
not be further considered, i.e., if 

Xj<X c f- a c - aj (67c) 
Yj < Ycf- a c - aj (67d) 
20 or 

Xj>Xd+a c +aj (67e) 
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Yj > Yd + a c + aj , (67f) 
then the j-th sphere lies beyond the roll corridor in the horizontal plane and is not a 
candidate for touching the current sphere. 
Detail of Subroutine FOBJ3 
5 Subroutine FOBJ3 finds the third object to be contacted by the descending 

current sphere when it is already in contact with two and only two other objects. 
There are only two possibilities for the two initial contacts: (1) both are placed pack 
spheres or (2) one is a placed pack sphere and the other is the cylinder wall for the 
submode of the current sphere. Each of these possibilities will now be addressed in 
10 turn. 

Touching two spheres 

It is presumed at the outset of this discussion that the two touching spheres, i 
and j, are in stable (compressive) contact, subroutine STABLE already having been 
called to remove them if they are not. A roll corridor follows a circle defined where 

15 the current sphere is in contact with both Spheres i and j. Of course, the current 
sphere cannot touch both spheres unless the condition rji < ai + aj + 2a c is satisfied 
where rji is the center-to-center distance between Spheres i and j. 

The current sphere's descent along this roll corridor will be defined by a 
single variable cp c '. This variable is the azimuthal angle of a local tilted coordinate 

20 system whose origin is in the center of Sphere i, with z axis tilted so it passes 
through the center of Sphere j, and with x axis swung upward so that it lies in a 
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vertical plane. In that coordinate system (which is denoted herein by primes), cp c ' 
alone describes the current sphere's position along the roll corridor. In the unlikely 
event that the current sphere finds itself perfectly balanced on Spheres i and j so 
that cpc' = 0, then a random number generator determines in which direction the 
5 current sphere will roll — toward negative (p c ; or toward positive <p c '. 

To find the coordinates of a pack sphere, say the k-th sphere, in this primed 
coordinate system, the preferred implementation first translates the global sphere 
positions to the unprimed local coordinate system whose origin is at the center of 
Sphere i but whose axes are parallel to the global axes: 
10 (68) 

Capital letters are used herein to denote coordinates in the global coordinate 
system, small letters in the unprimed local system (local to Sphere i), and primed 
small letters in the local tilted coordinate system. 

The Euler matrix used to rotate the unprimed into the primed coordinate 
15 system is given by 



A = 



^ - cos 0j. cos g>j. - cos 9 jt sin q> .. sin 9 }i ^ 
sin <pj. - cos (p n 0 



(69) 



sin 0ji cos <Pji sin 0 jt sin <p }t cos 9 }i j 

where 6ji and (pji are the polar and azimuthal angles of Sphere j relative to Sphere i 
and are given by 

cos#; = (Zj - Zi)/rji (70a) 
20 and 
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<pji = tani[(Yj-Yi)/(Xj-Xi)] . (70b) 
This Euler matrix allows the preferred implementation to obtain the k-ih sphere's 
position in the primed coordinate system as 

&' = A.Ji . (71) 

5 The current sphere's initial position in the primed coordinate system is given 

in like manner: 

£'=A-Ji . (73) 

Now that the three components of the vector I;/ have been obtained, the preferred 

10 implementation can find the initial azimuthal position 

qkd = tan^CydVoci') . (74) 



If cpci' is positive, the current sphere will roll toward higher positive values. If it is 
negative, it roll toward more negative values (higher in magnitude). In the unlikely 
15 case that it is zero, a random number generator chooses which direction it will roll. 

The polar angle of the center of the current spheres in the primed coordinate 
system is constant as long as it is in the roll corridor. It is obtained from the law of 
cosines as 

cos0c' = [rji 2 + (ac + ai) 2 - (a c + a;) 2 ] [2 rji (a c + a*)]' 1 (75) 
20 with 

sinGc' = +[1 - cos 2 6 c '] 1/2 . (76) 
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There are several possible outcomes as the current sphere descends its roll 
corridor. It can 

• hit the catch net or water level; 

• lose contact with Sphere i or Sphere j or both simultaneously; 
5 • hit the cylinder wall; or 

• hit another pack sphere or several pack spheres simultaneously. 

The roll corridor ends when the current sphere encounters the first of these 
situations. The preferred implementation then finds the q> c ' at which one of these 
events happens. The event with the smallest magnitude for q> c ' will be the bottom 
10 of the roll corridor and will be denoted cp c /. 

Hitting the catch net or water level 
First, the preferred implementation finds the global altitude of the current 
sphere as a function of cp c ': 



15 where the inverse Euler matrix is equal to its transpose A 1 = A T because of its 
orthogonality. The current sphere's global Z component is given by 



Its south pole hits the catch net if Zc = Solving for the primed azimuthal 

angle <p c r at which this equation may be satisfied, we have 



C I c 



(77) 



Z c = Zi + (etc + a;)(sin9ji sin0 c ' coscp c ' + cos0 3 i cosG c ') . 



(78) 



20 



cos = 




-- cos 0,, COS 0. ' 

Jl c 



(79) 



sin 0„ sin 0' 

J 1 c 
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where the sign of (pc is the same as the sign of <pd\ the angle at the top of the roll 
corridor. If this equation does not have a real solution in (0,7r/2), then the current 
sphere does not reach the catch net in this roll corridor. If <p c ' is outside the range 
[0,27r) ? then units of 2% are either added or subtracted from 9c' until it does lie 
5 within that range. Because neither Gji nor 6 C ' can ever be zero or n 9 Eq. (79) is never 
singular. 

The preferred implementation also can find which water level the current 
sphere might hit. Water level altitudes between adjacent cylinder walls are 
generally different as the pack is constructed. The only condition under which the 
10 current sphere might contact the k-th water level Z w (k) is if and only if the distance 
P c of the current sphere from the global Z axis lies between cylinder wall radii Wk-i 
and Wk at the point of contact. P c is given by 

Pc=(Xc 2 + Yc 2 y f2 (80) 
where X c and Yc can be obtained from Eq. (77) as 
15 X c = Xi + (a c +a;) (-cos 6fc cos^t sinft f coscpc + sin^i sinft' sin^b' 

+ sin$; cos^i cos ft') (81a) 

and 

Yc - Yi + (a c +ai)(-cos0ji sin^t sinft' cos<pc - cos#y; sinft' sin^fc' 

+ sin$; sin^i cos ft') (81b) 
20 Using the square of P c rather than Eq. (80), the water level checked for contact is 
Z w (k) if Wl_ x < Pi < Wl . Actually, each annular region between the (fe-l)-th and fe-th 
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cylinder walls is checked to see if the roll corridor has a possibility of contacting the 
water level there. 

The expression for the <p c ' at which the south pole of the current sphere 
touches the k-th water level is the same as Eq. (79) with Znet replaced by Z w (k). 
5 Then the end of the roll corridor <p c f is either ±n/2 or a contact point with the catch 
net or a water level as given by Eq. (79). 

Losing Contact With One or Both Spheres 
The preferred implementation also can determine the cp c ' at which the current 
sphere will lose contact with either Sphere i or j or both if it encounters nothing else 
10 as it descends its roll corridor. 

Of the two spheres, i and j, the one whose center is higher by the subscript u 
and the one whose center is lower is denoted by L Their global positions are (X u , Y u , 
Z u ) and (Xi, Yi, Zi) and their radii are a u and at, respectively. We always have that 
Z u > Zu Consider the unprimed coordinate system whose origin is in the center of 
15 the lower sphere and whose axes are parallel to the global axes. Let 6 C and cp c be the 
polar and azimuthal angular position of the current sphere in the unprimed local 
coordinate system centered in the lower sphere. 

As before, the primed coordinate system is the system obtained by tilting the 
unprimed z axis into the center of the upper sphere and then rotating the xy axes of 
20 that tilted coordinate system about its 2/ axis until its x'axis is pointing as steeply 
upward as possible, i.e., has a maximum positive projection onto the unprimed z 
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axis. This rotation is what Eq. (69) does if i is replaced by I and j by u. The polar 
and azimuthal angular positions of the current sphere in the tilted coordinate 
system are 0 C ' and q> c ', respectively where cos0 c ' is given by Eq. (75) with the same 
change in indices mentioned in the previous sentence. 



and upper spheres are at the same altitude, it will lose contact with them both 
simultaneously. To find the <p c ' position at which the current sphere loses contact 
with the upper sphere, the preferred implementation finds the two longitude lines 
of the lower sphere in the unprimed coordinate system along which the current 

10 sphere would just nick the upper sphere if it were descending along those longitude 
lines. The two <p c at which that happens are the ends of the roll corridor at which 
the current sphere would lose contact with the upper sphere and roll on the lower 
alone. Because the current sphere is rolling down only one side of the roll corridor, 
only that side's cp c is of interest. 

15 In the unprimed local coordinate system, the current sphere is at 

? c =(a c +a i )(ism8 c cos$ c + jsin# c sin^ c + kcos# c ) (82) 
and the upper sphere at 



5 



The current sphere will lose contact with the upper sphere first. If both lower 



r tt (isin# tt cos^ + jsin^ B sin^ B +kcos# w ) 



(83) 



Consider the equation 



20 




(84) 



If both sides are squared and written explicitly, we obtain 
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[(ac + ai) 2 sin9 c cosq) c - x u ] 2 + [(a € + at) 2 sin0 c sincpc - y«] 2 + 

[(a c + a/) 2 cos0c - z u ] 2 = (ac + a w ) 2 (85a) 

which can be rewritten 

r 2 +(a +a,) 2 - (a +a ) 2 

u 0{ \ \ c — = *u sin 9 C cos <p c + ^ sin 0 C sin ^ c + z u cos 0 C . (85b) 

5 The components x u , yu, and z u are given in Eq. (83). If they are substituted into Eq. 
(85b) and both sides are divided by r u , we obtain 

^, 2 +(g c +fl/) 2 -(flc +fl J 2 

2(fl c +fl,)r B (86) 
= sin 6> M cos <p u sin # c cos <p c + sin 0 W sin <p u sin 0 C sin #> c + cos 9 U cos 0 C 

where we recognize the left hand side as cos9 c ', the polar angle in the primed 
coordinate system which is constant as the current sphere descends the roll 
10 corridor. Now Eq. (86) can be written as 

cosGc' = sinGu sin0 c (coscp u coscp c + sincp u sincpc) + cosG u cos9 c (87a) 
which through a trigonometric identity can be written 

cosGc' = sinGu sinG c cos((p c - cp u ) + cosG u cosG c . (87b) 
We wish to solve for G c from this equation. For a given <p c , the G c solution has three 
15 possibilities: 

1. There are two Q c solutions if the current sphere hits the upper sphere as it 
rolls down the cp c longitude line of the lower sphere - one solution where it touches 
the upper sphere from above and the other where it touches from below; 
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2. There is one solution (the two solutions having become degenerate) where the 
current sphere just nicks the upper sphere while rolling along the <p c longitude line; 
and 

3. There are no real solutions, meaning the current sphere altogether misses 
5 the upper sphere while rolling along the <p c longitude line. 

The second possibility is the one of current interest because that is the point 
at which the current sphere can lose contact with the upper sphere and roll on the 
longitude line of the lower sphere. We will call that point the "departure point." 
Therefore, the preferred implementation first finds the angular position (0 C , cp c ) of 
10 the departure point and then finds the azimuthal angle q> c ' in the primed coordinate 
system to which this angular position corresponds. 

To solve for 0 C , given a particular <p c , the preferred implementation recasts 
Eq. (87b) is recast into the form 



cosfl c f 



= sin a sin 9 C + cos a cos G c 



(88) 




15 where 



sm a = 



sm0 u cos(<p c -<p u ) 



(89a) 




and 



cos 9, 



u 



(89b) 




Applying our trigonometric identity to the right hand side of Eq. (88), we obtain 
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10 



cos(0 c - a) = 



cos9 c } 



^sin 2 G u cos\<p c -<p u ) + cos 2 9 U 



(90) 



which gives the solution 



0 C = a + cos" 



cos & 1 



^/sin 2 0 W cos 2 (^ c - %) + cos 2 9 U 



(91) 



The arccosine generally has two distinct solutions: (1) the principal value which 
5 lies in the interval [0,7r] and which is denoted by capitalizing the function as Cos" 1 
and (2) the negative of the principal value. We are interested in the case where 
these two solutions are degenerate. That can only be so if the arccosine and its 
negative are the same value, i.e., zero. Therefore G c = a at the departure point, and 
the argument of the arccosine in Eq. (91) must be unity: 
cos ft ! 



^/sin 2 0 U cos 2 O c - q> u ) + cos 2 9 U 



■ = 1. 



(92) 



This equation is the condition used in the preferred implementation to find the 
unprimed azimuthal position of the departure point. Squaring both sides gives 
cos 2 e c ' = sin 2 9 u cos 2 ((p c - (p u ) + cos 2 e u (93) 
which can be solved for q> c to yield 



!5 <P c = <Pu + cos " 



/ cos^'-cos 7 ^" 
sin 2 6 U 



(94) 



The square root has a positive sign because the difference between cp c and cp u can be 
no greater in magnitude than nl2. Again, the arccosine has two solutions: the 
principal value Cos" 1 and its negative -Cos" 1 . The current sphere would lose 
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contact with the upper sphere if it is at azimuthal position <p u ± Cos' 1 . Its true 
position depends on where it was initially. If its initial position cp C i > 9u, then the 
principal value Cos -1 is used. If its initial position <p C i < <pu, then the negative of the 
principal value is used -Cos 1 . In the unlikely case that the initial position equals 
5 cp u , then the current sphere is balanced on the upper and lower sphere and a 
direction is randomly chosen by the random number generator. The preferred 
implementation using these relationships to determine where contact is lost. 

The preferred implementation also determines (0 C , cp c ) of the departure point. 
It then finds the departure point in the primed coordinate system by use of the 
10 Euler transformation of Eqs. (69) and (73). 

The first component of the primed vector is given by 
sin0 c ' coscpc' = - cos6 u sin0 c (coscp u cos<p c + sin<p u sincp c ) + sin9 u cos6 c (95) 
where the expression in parentheses can be written cos(cp c - <p u ). Eq. (94) can 
substitute for cos((p c - q>u). Also because 0 C is equal to a of Eqs. (89), we have 



Jsin 2 fl '-sin 2 ft 1 

15 sin 0 C = 1 u — e - (96a) 

cos 6^' 

cos# 

cos0 c = ^. (96b) 

cos 9 c y 

Substituting Eqs. (94) and (96) into (95), we obtain for the end of the roll corridor 
due to departure 

, cos#>n£ c ? tan^ ! 

ros «v Ws^-lift • <97) 
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The preferred implementation uses this relationship to obtain 0 C and q> c . Notice that 
the current sphere cannot be in stable contact with two pack spheres unless 9c < 6L 
Therefore q>cf always lies between 0 and 7t/2. Also, remember that cos9 c ' is given by 
the left hand side of Eq. (86) and that sinft' = +[1 - cos 2 &'] 1/2 . 
5 The second component of the primed vector is given by 

sin0 c f sincpc' - (sincp u coscp c - coscpu sin<p c ) sin9 c (98) 
and the same set of substitutions leads to 

J sin 2 6 » -sin 2 &' 

sin <p rf ' = . (99) 

^ cf cos6>'sin6> V } 

The 9c/ defined by Eqs. (97) and (99) is the departure point and is the end of the roll 
10 corridor. At this position, the current sphere loses contact with the upper sphere (or 
both if their centers are the same altitude). These relationships also are used in the 
preferred implementation to determine 0 C and <p c . 
Hitting the Cylinder Wall 
The cylinder wall for the current sphere is located at global radius W. The 
15 preferred implementation can find whether the roll corridor of the current sphere on 
the two pack spheres is such that the current sphere's center can intersect the 
cylinder wall. If it does, then the q> c ' of intersection is found. 

We need the global position of the current sphere as a function of <p c '. Again, 
the global position of the current sphere is 
20 1=%+% (100) 



93 



where J^. is the global position of one of the contact spheres and J; is the current 

sphere's position in the unprimed local coordinate system. (Remember the 
unprimed coordinate system is centered in Sphere i and has axes parallel to the 
global axes.) The unprimed vector is given in terms of the primed vector (that 
5 contains <p c ') by the inverse Euler matrix equation 

? =A- 1 -£ (101) 

where A~ 2 is the transpose of the Euler matrix in Eq. (69) and rc - (x c f , yc, Zc) is 
given by 

xc = (a c + ai) sinGc' coscp c ' (102a) 
10 y c r = (a c + at) sinG c ' sincp c ' (102b) 
zc 1 = (a c +ai) cosGc' . (102c) 
The preferred implementation therefore finds the global position R c = (Xc, Y c , Z c ) in 
terms of the single variable cp c ' describing the descent down the roll corridor. If the 
lateral global position of the current sphere reaches the cylinder wall radius Wfox a 
15 value of (p c f in the interval (<p C i', q>c/), then the current sphere can reach the cylinder 
wall if it does not encounter another object first. In other words, if there is a 
solution of 

Xc 2 + Yc 2 = W (103) 
for q> c ' in the interval (cp C i', q>cfO> then the current sphere reaches the wall. Written 
20 explicitly, Eq. (103) has the form 

Ao +Ai coscpc' + As cos 2 cpc + Bi sincp c ' + B2 sin 2 cp c ' = 0 (104) 
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where 

Ao = (Xi 2 + Yi 2 - W)l(a c + at) + sinGjt cos9 c ' [2Xi coscp;; + 

2Yi sincpyj + (a c + ai) sinG/i cosG c '] (105a) 
Ai = -2 sin0 c ' cosGyt [Xi cos(p>i + Yi sincpjj + (a c + ai) sinG;i cos9 c '] (105b) 
5 A2- (a c + at) cos 2 9y* sin 2 9 c ' (105c) 
Bi ~ 2 sinGc' [Xi sincpyi - Yi coscpyj] (105d) 
B 2 = (a c + a») sin 2 G c ' . (105e) 
Eq. (104) can be converted into a quartic equation in sin(p c ' by changing the cos 2 cp c f 
term into 1 - sin 2 <p c ' and the coscpc' term into +[1 - sin 2 (p c '] 1/2 - (The plus sign only is 
10 used only because <p c ' cannot exceed ±tc/2.) We then have 

(Ao + A 2 ) + Bi sinq>c' + (Bz - Ai) sinV = -Ai[l - sin 2 <p c '] 1/2 . (106) 
Squaring both sides gives 

(As - B2) 2 sinV - 2Bi(A 2 ~ B2) sinV + (Ai 2 - 2A0A2 - 2A2 2 + Bi* + 2AoB 2 + 2A2B2) 
sin 2 ^' + 2(Ao +A 2 )Bi sin^ c ' + (Ao - Ai + A 2 )(Ao + A1+A2) =0 

15 (107) 

If this equation has any real roots in the interval ((p C i\ q>cfO> then the smallest in 
magnitude of these that satisfies the original equation (104) is the point at which 
the roll corridor intersects the cylinder wall. The preferred implementation can 
check the roots in the original equation because squaring can give extraneous roots. 

20 Hitting One or More Spheres Simultaneously 
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The preferred implementation has now defined a roll corridor in the primed 
coordinate system beginning at <p C i and ending at qkf. The end of the roll corridor 
is where the current sphere either lost contact with one or both of its contact 
spheres, when its south pole encountered either a catch net or a water level, or 
5 when its center encountered its cylinder wall. The preferred implementation can 
check all spheres that can intersect this roll corridor to find the sphere(s) 
encountered first as the current sphere descends along its roll corridor, i.e., as <p c f 
moves from <p& to <p c f. 



10 (77, 78, and 81). If some other pack spheres, say the fe-th sphere is to intersect the 
roll corridor, then the following equation must have a solution in the range {q>ci,(pc{)\ 



The global position of the current sphere as a function of <pc is given by Eqs. 




(108) 



Squaring this equation gives 



R 2 c -2$ c l k +Rl={a c +a k ) 2 



(109) 



15 which becomes the following equation in <p c f : 



A cos^' + B sm<pc = C 



(HO) 



where 



A = -X ik cos <p~ cos 0j. sin 6 C '-7 ft cos 9 fi sin <p }i sin 9 C '+Z ik sin 0 C ' sin 6 }i 



(HI) 



B = X ik sin <p fi sin 9 C '-Y ik cos <p Jt sin 9 C ' 



(112) 
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C--X ik cos pj. cos 0 C ! sin 0 yi -Y ik zo%6 c y sin <p jt sin 0.. - Z ik cos # c ! cos 
(2a c +g f +aj(^-a,)-^ 



(113) 



+ 



2(a c +fl,) 

The global coordinates of Sphere & relative to Sphere i are given by 
Xik = Xi-X k , Yik = Yi-Yk, Zik = Zi-Zk, (114a) 
and the square of the distance between the two sphere centers is 
5 Rik 2 = Xik 2 + Yik 2 + Zik 2 . (1 14b) 

The solution procedure of Eq. (110) is similar to that given in Section VILD.7. The 
preferred implementation finds (pc to be 

^ =a±C0S ''7T^W (115a) 

where 

10 a = tan -1 -. (115b) 

If any of the <pc r solutions are outside the range [0,27i), units of 2tc are added or 
subtracted until they are inside that range. Then if any solutions lie within 
(<pti',q>cf), the smallest in magnitude of these is the position of first contact between 
the current sphere and Sphere k. If there is no solution, then Sphere k intersects no 
15 part of the roll corridor. 

The preferred implementation has now performed an exact solution to find 
whether and where Sphere k intersects the current sphere's roll corridor. In most 
packs, the vast majority of the pack spheres are out of range of the roll corridor and 
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so do not need the detailed calculations of Eqs. (115). There are a few much simpler 
checks that can be performed by the preferred implementation to immediately 
determine whether a pack sphere is out of range and so is not a contact candidate. 
First, the roll corridor's highest global altitude is at qki which, from Eq. (78) 
5 is found to be 

Zc,max = Zi + (a c + aO(sin0ji sinG c ' coscpci' + cosfyi cos0 c ') - (116) 
Therefore, the k-th pack sphere whose center lies a distance a c + ak above Z c ,max 
cannot overlap the roll corridor. Likewise, the lowest global altitude is at <p c f and is 
given by 

10 Zc 9 min = Zi + (a c + a0(sin6ji sin0 c ' cosepe/ + cosOji cos9 c ') . (117) 

The k-th pack sphere whose center lies a distance a c + ak below Zemin also cannot 
overlap the roll corridor. Therefore, the preferred implementation can immediately 
exclude from consideration all spheres whose global altitudes Zk satisfy the 
following: 

15 k-th sphere does not overlap if 

Zk > Zc,max + a c +ak (118a) 
or 

Zk < Zc } min - a c - ak . (118b) 

Roll corridor extrema in the Xand Y directions do not necessarily lie at the 
20 beginning and end of the roll corridor. Eqs. (81) give the global Zand Y coordinates 
of the current sphere as a function of If we differentiate X c (<pc) with respect to 
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q>c and set the derivative equal to zero, we have an equation whose solution gives us 
the extrema in X. If either of those extrema lie in (<pti\(pcf), then that position is 
either the minimum or maximum X value the current sphere assumes as it 
descends the roll corridor. The same is true for the derivative of Y(qk f ) with respect 
5 to <p c \ The derivatives followed by the extrema solutions are 



d<p t 



7 = (a c + a t )(cos 0 }i cos <p fi sin 0 C ' sin <p c 4- sin <p jt sin 0 C ' cos <p c ' ) = 0 (119) 



tan <p {; 

tan^*=-— f- (120) 

COS &j . 

for the two X extrema (the above equation has two solutions for qhx% and 
= (a e + a t )(cos 0j. sin (p }i sin 0 C ? sin (p c - cos ^ sin # c ' cos (p c ! ) = 0 (121) 

/C 

cote?,, 

10 tan^'=— ^ (122) 
' cosfcL 

for the two Y extrema. 

Finally then, the X and Y extrema of the roll corridor in the global 
coordinates are 

Xc,min = MIN[Xc(^/), Xc(<Pcf), Xc(<pcx')] (123a) 

15 Xc,max = MAX[X C (^'), M<Pcf), M<P*')] (123b) 

with the Xc(<pcx) omitted if neither (pal lies in {<pd' ,<pcf). Similarly, 

Yc,min = MIN[Yc(^£'), Y c (<&/), Y c (qky')] (124a) 

Yc.max = MAX[ Yc(pbi'), Y c (<pcf), Y c (^')], (124b) 



99 



again with the Y c {(pcy) omitted if neither q> C y lies in (<pri f ,<pef). 

Then the fe-th sphere cannot overlap the roll corridor if 
Xk < X c ,min - Oc- ak (125a) 
or 

5 Xk > Xc,max + a c + ak (125b) 
or 

Yk < Yc,min - a c - ak (125c) 
or 

Yk > Yc,max + dc+Clk. (125d) 

10 Another discriminator that can be used by the preferred implementation to 

immediately exclude a large number of pack spheres from detailed examination is 
to consider whether the roll corridor is far enough from the central Z axis that it 
cannot overlap all the spheres within a given inner cylinder. To find this, the 
preferred implementation uses the closest radial approach P c ,min of the roll corridor 

15 to the central axis. The square of the radial distance [Pc(<pe')] 2 is given by 

[Pc(<p c 'W = [M<pc'W + [Yc((pc r )] 2 (126) 
where we obtain the explicit expression by substituting Eqs. (81) into Eq. (126): 
[Pc(<pc)] 2 = [Xi + (cic^Gi) (—cos 0ji costpji sin&' cos<pc f + sin^ sin&' sin^' 

+ sin&ji COS^r COS&')] 2 

20 + [Yi + (ctc+CLi) (—cos Oji sin0$i sinft' cos<p c f - cos^'i sinffc sin^c' 

+ sin#; smyji cos&')] 2 . (127) 
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Setting the derivative of [Pc(<pc)] 2 with respect to cpc equal to zero yields the 
equation whose solution <p C p gives the closest approach of the roll corridor to the 
global Z axis: 

[Xi + (ac+cti) (-cos fyi cos^j sinft' cosg?e + sin^i sinft' sin^ f + 
5 sin^fj cos#u cos ft')] (cos 0ji cos#$; sinft' siaqk' + sin^i sinft' cos<p c r ) + 
[Yi + (a c +ai)(-cos0ji sin<pji sinft' cos<pc - cos^i sinft' sin<pc r + 
sin^i sin^t cos ft')] (cos 0ji sin^i sin 9c sin^c' - cos^j sinft' cos#> c ') = 0 . 



(128) 



This equation is of the form 



10 Ao +Ai cos<p c r + A2 cos 2 <p c f + Bi sin^' + B2 sin 2 ^' = 0 



(129) 



where in this case, the coefficients are given by 



Ao - Xi 2 + Yi 2 + 2(a c + ai) cos ft' (X sin#; cos#/; + Yi sinlfc sin^i) + 



(a c + ai) 2 cos 2 ft' sin 2 ^fj 



(130a) 



Ai = -2(a c + ai) cos0ji sinft' [Xi cos<pji + Yi sin<pji + (ac + ai) cosft' sin$i] 



15 



(130b) 



A2 = {a c + ai) 2 cos 2 &ji sin 2 ft' 



(130c) 



Bi = 2(a c + ai) sinft' {sin 6$; [Xi + (a c + ai) cosft' sin#; cos^i] 



cosft'[Yi + (a c + ai) cosft' sinlfc sin^yi]} 



(130d) 



20 S 5 = (ac + ai) 2 sin 2 ft' 



(130e) 
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Eq. (129) can be converted into a quartic equation by changing the cosines into 
sines in the same fashion as led up to Eq. (107): 

(As - B2) 2 sinV - 2Bi(A2- B2) sinV + (Ai% - 2A0A2 - 2A2 2 + B1 2 + 2A0B2 + 2A2B2) 
sinV + 2(Ao +A 2 )Bi sin^' + (Ao - Ai + A2XA0 + A1+A2) - 0 . 
5 (131) 

There are four sintpc solutions to this quartic equation and for a given sin<p c f , there 
are two q>c solutions. Complex solutions and real solutions lying outside [-1,1] are 
immediately discarded in the preferred implementation. Furthermore, we are only 
interested in those solutions with the same sign as %v&qki and that lie in ((pci,(pc{) 

10 because that is the direction the current sphere will roll. The preferred 

implementation checks any remaining solution <p C p by ensuring it satisfies the 
original equation (129). If there is no solution in [<pci r ,<pcf], then one or the other of 
the end points is the closest approach to the global Z axis and both must be 
substituted into Eq. (127) to find the smaller P c . 

15 Once the P c ,min of closest approach is found, the preferred implementation 

checks to see if there are any mode's cylinder radii Wk that lie completely inside 
Pcmin by more than the sum of the mode sphere's radius plus the current sphere's 
radius. If 

P c > Wk+ai + ac (132) 
20 then none of the Mode k particles can touch the current sphere as it descends its 

current roll corridor and so may immediately be excluded as candidates for contact. 
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All spheres whose north poles have sunk beneath the catch net or water level 
are not candidates for contact with the current sphere. The preferred 
implementation continually or periodically monitors the lowest index number of the 
sphere whose north pole is still above the catch net so that we waste no time 
5 looking at spheres with lower index numbers as candidates for contact. 



Touching One Sphere and the Cylinder Wall 

Remember that in the preferred implementation it is the center of the 
current sphere constrained by the cylinder wall, not its edge. When in stable 

10 contact with both a pack sphere and the cylinder wall, the roll corridor lies on the 
wall and descends toward the equatorial plane of the contacted pack sphere. If the 
contacted sphere (denote it Sphere i) is as large as the cylinder, then the current 
sphere moves toward the lowest point on a roll corridor that spans the entire 2% 
radians of the cylinder. 

15 The roll corridor will be denoted by the global azimuthal position <P C . Its 

initial position is given by 

0co = tan^Yc/Xc) . (133) 
The preferred implementation uses the range [0,2n) rather than (-n 7 7z\ for the range 
of the global azimuthal position. The global azimuthal position of the Sphere i is 
20 0i = tem'^Yi/Xi) (134) 

which also lies in [0,2n). If @ c o > 01, the current sphere rolls toward higher 0. If 0 c o 
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< 0i, it rolls toward lower 0. In the unlikely event that 0 c o = 0i, the random 
number generator chooses a roll direction. 

If it encounters nothing else, the current sphere will roll until it reaches: 

1. the equatorial plane of Sphere i where it will fall away from Sphere i; or 

2. the lowest azimuthal point in its cylinder permitted by Sphere L 

The azimuthal position at which it will reach the equatorial plane of Sphere i is 
found by the law of cosines as 

0cf= 0i±A0 c (135) 



In this expression, W c is the radius of the cylinder wall, P. = ^jxf + is the radial 

position of Sphere i from the global Z axis, and a c and ai are the radii of the current 
sphere and Sphere i, respectively. The plus sign is used in Eq. (135) if the current 
sphere is rolling toward increasing 0, the minus sign if toward decreasing 0. If the 
absolute value of the quantity in the parentheses of Eq. (136) is greater than unity, 
then there is no solution and the current sphere simply rolls around the perimeter 
of the cylinder until it reaches the lowest point it can. This lowest point is opposite 
the azimuthal position of Sphere i. Hence if 



where 




(136) 




W?+P>-(a e +a t ) 2 



(137) 



2W c P t 
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then 

A&c = ±tc . (138) 
In the unlikely event that Pi = 0, i.e., Sphere i is centered on the global Z axis, then 
there is no downward direction. The roll corridor is horizontal and the current 
5 sphere is placed where it is. 

In this analysis, if the addition or subtraction of A0 C causes 0 c fto lie outside 
[0,2tt), then units of 2n are added or subtracted until 0 C f again lies in [0,2tt). 

The preferred implementation also determines what placed pack spheres 
might intersect the roll corridor in the range (0 c o, @cf). 
10 The current sphere will impact a placed pack sphere, say the j-th sphere, if 

there is a solution to 

|^(O c )-^| = fl c +a y (139) 

anywhere in the range (0 c o,0cf). The Cartesian components of the vector 
K e (O e ) = (X C J C ,Z C ) are given by 

15 Xc(0c) = Wc COS0c (140) 

Y c (0c)= W c sin0c (141) 
The Z c component can be found by inserting the X c and Y c component into the 
condition 

|£-J(| = a c +a,. (142) 

20 Squaring both sides, we obtain 

Wc 2 + Zc 2 + Ri 2 - 2W c (Xi co$0c + Yi sin0 c ) - 2ZiZ c = (a c + at) 2 (143) 
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which is quadratic in Zc. One of the solutions is when the current sphere is above 
Sphere i and the other when it is below. In the preferred implementation, we are 
only interested in the solution above Sphere i, which is 
Zc(0c) = Zi+ [Zi 2 + (etc + a;) 2 - W c 2 - R? + 2W c (Xi cos0 c + Yi sin<Z> c )]i/2 . (144) 
5 Now, inserting Eqs (140, 141, 144), in Eq. (139), we obtain an equation to solve for 
the 0c at which the current sphere contacts the j-th pack sphere: 
Ao + Ai cos0c + Bi sin 0c + (A2 cos0 c + B2 sin<Z> c ) 2 = 0 (145) 
where the coefficients are given by 

Ao = G 2 - 4(Zi - Zj)2[Z? - Wc 2 - ifc2 + ( ac + ai )2] ( 146a ) 
10 Ai = AGWcQCi - Xj) - 8(Zi - Z^WcXi (146b) 
A 2 = 2Wc(Xi-Xj) (146c) 
Bi = 4GW c (Yi - Yj) - 8(Zi - Zj) 2 W c Yi (146d) 
B 2 = 2W c (Yi-Yj) (146e) 
and G is given by 

15 G = (m - aj)(2a c + ai + aj) - Ri 2 + Rj 2 + 2Z? - 2ZiZ, . (147) 
Eq. (145) can be converted into a quartic polynomial in cos<Z> c : 

(As 2 + B2 2 ) 2 COS 4 0c + 2(AlA2 2 + 2A2BlB2 - A1B2 2 ) COS 3 0c + 

(Ai 2 + 2A0A2 2 + Bi 2 - 2A0B2 2 - 2Ab 2 B2 2 - 2B2 4 ) cos 2 0 C + (148) 
2(AoAi - 2A2B1B2 + A1B2 2 ) cos 0c + (Ao -Bi+ B^Ao +Bi+ B2 2 ) = 0 
20 whose solutions are candidates for contact of the j-th pack sphere. The preferred 
implementation uses this expression as part of its determination. There are four 
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solutions of the quartic equation and each cos0 c solution has two 0c solutions — 
eight in all. Of the four quartic solutions, the preferred implementation can 
eliminate complex solutions or solutions whose magnitude is greater than unity. Of 
the remaining solutions, if there are any that satisfy the original equation (145) and 
5 that lie in the range of the roll corridor (@ c o,&cf), they are contact points. Of the 
contact point(s), the value of 0 C that lies closest to 0 c o is the point at which the 
current sphere touches the j-th sphere. We will denote that point by @ C j. If there 
are no solutions that satisfy the above conditions, the current sphere cannot contact 
the j-th sphere. 

10 As the preferred implementation checks the candidate spheres in the pack for 

contact with the current sphere, that sphere with 0 c j closest to 0 c o, if any, is the 
next sphere contacted. In the unlikely case that two or more pack spheres have the 
same 0 c j, then it checks the stability of each trio of contacts to see if the current 
sphere is to be placed there. If it is not stable, then the stability routine determines 

15 which sphere(s) the current sphere will roll away from. 

As an efficiency measure, the preferred implementation can eliminate from 
consideration any modes or categories that are sufficiently far inside the cylinder 
wall of the current sphere that they cannot touch it. All mode k spheres are out of 
range if 

20 Wk + Gk£Wc-ac. (149) 
Any other sphere j is a candidate for being in range only if 
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Wc-ac- aj <Pj<Wc + a c +aj (150) 
where Pj is the global radial distance of Sphere j from the Z axis. 



THE SUBROUTINE UPDATE 

5 The UPDATE subroutine of the preferred implementation updates all 

parameters and arrays that change as a result of the placement of the current 
particle. These changes include the following. 

1. The Pool Switch 

If the pool switch is true, then the placement of the current particle means 
10 the bin from which it was chosen should be diminished by one. This adjustment is 
performed in Subroutine UPDATE. 

2. The Catch Net 

The catch net belonging to the mode, submode or category current sphere also 
generally ratchets upward. The prescription for determining how much it ascends 
15 will now be addressed. 

A sense of pack surface must be established. With a variety of particles being 
dropped into a variety of cylinders, the pack surface must be defined. This can be 
done in a number of ways, as noted above. For purposes of the preferred 
implementation, let us consider the altitude of the south poles of the last m placed 
20 spheres of any and all modes. As explained above, the average of those altitudes is 
defined for present purposes as the pack "surface": 
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1 



m 



'surface 



m 



VILE.2(1) 



where m is the number of most recently placed pack sphere equatorial cross- 
sections that would just cover the outermost cylinder n with cross-sectional area 
nWn 2 . The value for m is found from the inequality 



If the size of the last placed set of spheres is very large so that m is very small, then 
we have a minimum required value for m which can be user-defined but for which 
the default value in this implementation IS TYlmin = 12. Only when the pack is just 
starting so that 12 particles have not yet been dropped is m smaller for this 
implementation. Pre-placed particles do not count in assessing the pack "surface." 
The &-th mode catch net is then defined as 

Znet(k) = Zsurface ~ jak VTLE.2(3) 

where y is the number of £-th mode radii beneath the surface that the £-th mode 
catch net lies. Here it is a user-defined number which has a default of y - 2. 

Subroutine UPDATE readjusts all the above parameters so that the catch net 
for every mode has a proper value after the current sphere is placed. 

When the pack is still small enough that an insufficient number of particles 
have been placed to satisfactorily define a surface, the catch net may be ratcheted 
up a little after the placement of each particle so that the particles on the bottom of 



m-\ m 



VII.E.2(2) 
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the cylinders do not all have the same altitude. If the switch ROUGHNET is true, 
then the current sphere's catch net moves upward by a small random amount 
ramiJNcs VILE.2(4) 
where r is a uniformly generated random number in [0,1] and Ncs is the total 
5 number of sphere cross sections required to cover all the cylinders: 
« w 2 - w 2 

2 Vn.E.2(5) 

k=\ a k 

where n is the number of modes, Wk is the cylinder wall of the fe-th mode, and Wo = 
0. If ROUGHNET is false, then the catch net stays where it is until enough 
particles have been defined to make a surface and until the surface is more than 

10 Zsurface — JQ>k above the starting place. 

The starting place for the k-th mode catch net is Zinit - ak. A fe-th sphere 
radius is subtracted from the overall pack starting altitude Zinit because the catch 
net stops south poles and Zinit is the starting place for sphere centers. The default 
for Zinit is zero in this implementation. 

15 The Water Levels 

The water level also is modified when the current sphere is placed according 
to the preferred implementation. Recall that the water level only exists beyond the 
innermost cylinder. Its altitude tracks the pack surface and is at a user-specified 
distance p in terms of the smallest sphere radius amin beneath the surface defined 

20 by Eqs. VILE.2(1) and (2). The default for p is unity. Hence 

Zw — Zsurface ~~ Pdmin . VILE.3(1) 
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If p is negative, the water level lies above the pack surface. 

The smallest cylinder in this implementation has no water level, only a catch 
net. It does not need a water level because water levels are designed to simulate 
the surface of the innermost cylinder and the innermost cylinder has its own 
5 particle surface potentially containing all the particles of the pack. 
Current Sphere History Parameters 
The preferred implementation uses are parameters that monitor the number 
of times the current sphere tunneled, the number of distinct roll corridors it 
underwent (not counting free-fall drops), and the number of free-fall drops it 
10 underwent. All these are reset to zero in this implementation. 

The arrays indicating which pack spheres the current sphere is touching and 
which spheres it just rolled away from are also zeroed out along with the total 
number of such spheres. 

The array keeping track of the lowest index number above catch net for each 
15 mode is adjusted if necessary. 

The array monitoring the index numbers of all pack spheres with exposed 
north poles also is updated if the current sphere covered one or more exposed north 
poles. 

In addition, the pack sphere position arrays have the current sphere's (X,Y 9 Z) 
20 position added to them along with the array specifying its mode number. The 
counter monitoring the total number of spheres placed so far in the pack is also 
incremented by one. 

Ill 



THE SUBROUTINE OUTSTAT 

There are a number of important pack statistics and features that the user 
may need or want to calculate. In this implementation, subroutine OUTSTAT 
calculates the output statistics that the user specifies. These output statistics 
include: 

A. The volume fraction of the pack. 

B. The radial distribution functions gij(r)dr which is the mean number of Mode j 
particle centers that lie between r and r + dr from a Mode i particle center. 

C. Coordination numbers ft, which represent the mean number of Mode j 
particles touching a Mode i particle. 

D. Detailed pack statistics for each particle mode such as that mode's minimum 
and maximum positions, number of its particles in the pack, radial and axial 
density functions for each requested mode, as well as other information. 
Each of these output data will now be addressed. 

Finding the Volume Fraction 

To determine the packing fraction, the program finds the total volume of the 
innermost cylinder occupied by spheres. Recall that the cylinder walls are 
boundaries to the centers of the spheres, not their edges. Hence, many spheres will 
be overlapping the cylinder walls. The preferred implementation has a general 
expression for determining how much volume of a sphere of radius a s located at 
position (x s , y s> z s ) lies inside a cylinder which is concentric to the z axis with radius 
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etc, bottom cylinder altitude Zb and top cylinder altitude zu The sphere volume 
inside the cylinder will be denoted Vin. Here the sphere and cylinder geometric 
parameters have no restrictions except: 

a s > 0, a c > 0, and zt > zb . VIILA.(l) 
5 The radial distance of the sphere center from the z axis is defined here by 

ps=[xs 2 + ys 2 ] m . VIILA.(2) 

There are four possible cases for this overlap problem. In Case 1, the sphere 

is completely outside the cylinder so Vin - 0. In Case 2, the sphere is completely 

inside the cylinder so Vin = (4/3)7ia s 3 . In Case 3, the cylinder is completely inside the 
^ 10 sphere so Vin = na c 2 (zt - zb). In Case 4, the sphere and cylinder intersect each other 
~?t and the expression for Vin is more complex. In the preferred implementation, 
q processing is performed to identify and address these cases using the following 
s approach. 

p; 1 Case 1: Sphere Completely Outside Cylinder 

f 15 The sphere is completely outside the cylinder 

Vn = 0 VIILA.1(1) 

only if one of the following conditions are met: 

If the sphere's north pole is lower than the cylinder bottom . . . 

z s + as<zb VIII.A.l(la) 
20 Or if the sphere's south pole is above the cylinder top . . . 

Zs — <Zs> Zt VIII.A.l(lb) 
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Or if the sphere is more than its radius away from the cylinder side . . . 
ps-a s >ac VIILA.l(lc) 
Or if the sphere center is above the cylinder top and is more than its radius away 
from the upper cylinder edge. . . 
5 (ps - ac) 2 + (z s - zi) 2 > as 2 andz s >zt VIII.A.l(ld) 
Or if the sphere center is below the cylinder bottom and is more than its radius 
away from the lower cylinder edge . . . 

{ps - a c ) 2 + (zb - zs) 2 > as 2 and z s <zb. VIILA.l(le) 

ilO Case 2: Sphere Completely Inside Cylinder 

l The sphere is completely inside the cylinder 

j Vin = (4/3W VIILA.2(1) 
only if all three of the following conditions are simultaneously met: 

I If the sphere's north pole is no higher than the cylinder top . . . 

jl5 zs + a s <zt VIII.A.2(la) 
And if the sphere's south pole is no lower than the cylinder bottom . . . 
zs-a s >z b VIII.A.2(lb) 
And if the sphere's equator nowhere exceeds the cylinder radius . . . 
A + a s <a c . VIII.A.2(lc) 

20 Case 3: Cylinder Completely Inside Sphere 

The cylinder is completely inside the sphere 
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Vin = nacHz t -z b ) VIII.A.3(1) 

if both the following two conditions are met: 

If the upper edge of the cylinder is inside the sphere . . . 

{ps + a c ) 2 + (zt - zs) 2 < as 2 VIILA.3(la) 

And if the lower edge of the cylinder is inside the sphere . . . 

{ps + ac) 2 + (zs - z b ) 2 < as 2 . VIILA.3(lb) 



Case 4: Cylinder and Sphere Intersect 

There are several subsets of Case 4. To efficiently organize them, we define 
two circles: the sphere circle and the cylinder circle. The sphere circle is the 
perimeter of the sphere's cross section at a given z altitude. Its radius is 
r s = [as 2 -(z- z s ) 2 ] 1/2 , zs - a s < z < z s + a s VIII.A.4(1) 
and its radial position ps is given by Eq. VIII.A(2). Likewise the cylinder circle is 
the perimeter of the cylinder's cross section which has constant radius a c at any z 
altitude between zb and zt and is centered on the z axis. 

To describe the intersection of sphere and cylinder, four subcases for the 
relationship of the sphere circle and the cylinder circle at a particular altitude z are 
defined: 

Subcase 4.1: The two circles do not overlap or one or both are nonexistent. 
Subcase 4.2: The sphere circle is completely enclosed in the cylinder circle. 
Subcase 4.3: The cylinder circle is completely enclosed in the sphere circle. 
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Subcase 4.4: The two circles intersect. 
To find Vin, the preferred implementation integrates in z the overlap areas of the 
two circles Ain from the altitude where one subcase holds to the altitude where 
another subcase holds. 
5 For Subcase 4.1, Ain = 0. VIILA.4(2a) 
For Subcase 4.2, Ain = na s 2 . VIII.A.4(2b) 
For Subcase 4.3, Ain = nctc 2 . VIILA.4(2c) 
For Subcase 4.4, J±in — Ctc 2 {a - cosa sina) + r s 2 (fi - cos/? sin/?) VIII.A.4(2d) 
where r s is given by Eq. VTILA.4(1) and 

2,2 2 

10 cosa = Ps c ~ r * VIII.A.4(3a) 

2 Ps a c 



Vfc 2 ~(« c -P s ) 2 ][(a c +A) 2 ~rf] 

sin a = Vill.A.4(3b) 

2 2 2 

cos6= Ps +r * ~ Uc VIII.A.4(4a) 
IPS, 



MlEz^mi^zSl . VIII.A.4(4b) 

The implementation now constructs a series of altitudes at which the 
15 cylinder circle and the sphere circle change from one subcase to another in Case 4. 
A sphere and a cylinder side can intersect each other at either 0, 2, or 4 
altitudes. Pairs of altitudes may be degenerate (two solutions having the same 
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altitude) if the sphere just nicks the cylinder wall. The preferred implementation 

determines these altitudes using the following equation: 

Mz)] 2 = [oc ± /*] 2 VIILA.4(5a) 

or 

5 a s 2 -{z- zs) 2 = [a c ± ps] 2 VIILA.4(5b) 
which, using Eq. (9), has the four solutions 

z = z s ± y s -(a c ±p s f . VIII.A.4(6) 

The + sign in the first ± of Eq. VIILA.4(6) denotes the upper or northern hemisphere 
intersection altitude; the - sign denotes the lower or southern hemisphere 

10 intersection altitude. The + sign in the second ± of Eq. VIII.A.4(6) denotes 

intersection with the far cylinder wall, on the opposite side of the z axis from the 
sphere's side; the - sign denotes intersection with the near cylinder wall, on the 
same side of the z axis as the sphere. 

Armed with these four solutions, the implementation can effectively address 

15 each of the four subcases of Case 4. In this illustrative implementation, it 

separately addresses whether the sphere circle lies in the southern hemisphere of 
the sphere or whether it is the equatorial plane or in the northern hemisphere. The 
four subcases are considered separately for the two hemispheres making eight 
situations in all. The hemispheres are considered separately because if the sphere 

20 circle is in the southern hemisphere, it is increasing in size as z increases. Hence 
the altitude at which we will change subcases can be determined. If the sphere 
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circle is the equatorial plane or in the northern hemisphere, it is decreasing in size 
as z increases and again the altitude at which subcases will change can be 
determined. 

Finding Upper and Lower Boundaries of Subcases for Case 4 



subcases change will now be discribed. Let ziow be the lower altitude of the subcase 
currently being examined and zugh be its upper altitude. The function r s (z) is given 
in Eq. VIILA.4(1). 

Initially, let ziow = zb, the bottom of the cylinder. After each DO WHILE cycle, 
10 ziow is set to the previous cycle's zhigh and a new Zhigh is sought. To illustrate using 
pseudocode expression: 

Vin = 0 ! Initialize Vin to zero. 

ziow = Zb ! Initialize ziow to the bottom of the cylinder. 

keepGoing = TRUE ! Do Loop flag. 

15 DO WHILE (keepGoing) 



5 



The logic used in the preferred implementation to find the altitudes at which 



IF (z s > ziow) THEN 



! Southern hemisphere subcases 



IF O - r s (ziow) > ac) THEN 



! Subcase 4.1 




IF (zhigh > zt) keepGoing = FALSE 



20 



(There is no addition to Vin in this subcase.) 



ELSE IF {ps + r s (ziow) < a c ) THEN 



! Subcase 4.2 
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zhigh = MIN[z s - ^a] ~(a c - p s f , zt, z s + a s ] 

Vin = Vin + n(zhigh - Zlow)[a,s 2 - Zs 2 + Z s (zhigh + £Zou>) - 
( Z AigA + 2 high 2 low + Z L)^] 

IF (^tgA = OR = ^ + as) keepGoing = FALSE 
ELSE IF (rs(ziow) - > ac) THEN ! Subcase 4.3 



zttgft = MIN[z, +<Jal~(a c +p s ) 2 , zt] 

Vin = Vin + nac 2 (Zhigh - Zlow) 

IF (sa#a = 2&) keepGoing = FALSE 
ELSE IF ((r s (ziow) > ps AND r s {zhw) - fk<a c AND ac < n{ziow) + /*) 
10 OR 0* > r s (ziow) AND ps - r s (ziow) < a c AND a c < ps + r a (zu>w))) 

THEN ! Subcase 4.4 

IF (as > a c + a) THEN 



zugh = MIN[z, ~^a 2 s -(a c +p s f , a] 

^ = P*. + = Vi» + AVi VIILA.4(7) 



15 where LVi is given in the next section below. 

IF (zugh - zt) keepGoing = FALSE 

ELSE 



zhigh = MIN[z s +^a 2 s -{a c -p s f , zt] 
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V in = V in + \ dzA in = Vin + AV 2 VIII.A.4(8) 

where AV2 is given in the next section below. 
END IF ! End two Subcase 4.4 situations 

END IF ! End Subcase 4.4 

ELSE ! Equatorial plane & northern hemisphere subcases 

IF (ps - r s (ziow) > ae) THEN ! Subcase 4.1 

keepGoing = FALSE 

(Because the sphere circle is shrinking away from the 
cylinder circle, there is no chance for additional overlap.) 

ELSE IF (p s + rsiziow) < a e ) THEN ! Subcase 4.2 

Zhigh = MIN[2f, Zs + a s ] 

Vin = Vin + n(zhigh - Zlow)[a s 2 - Zs 2 + Zs(zhigh + Slow) - 

( z ii g h + z high z tow + z h w y3] 

keepGoing = FALSE 
ELSE IF (r s (ziow) -ps> a c ) THEN ! Subcase 4.3 

zhigh = MIN[z s + 4 a2 s~i a c + Psf > Zt ] 

Vin = Vin + %ac 2 (Zhigh - Zlow) 

IF (zkigk = zi) keepGoing = FALSE 
ELSE IF ((r s (ziow) > p s AND r s (ziow) - ps<a c AND a c < r s (ziow) + ps) 
OR (p s > r s (ziow) AND p& - rsiziow) < a c AND ac < ps + r s (ziow))) 
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THEN 



! Subcase 4.4 



::10 



zhigh = MW[z s + ^a 2 s -{a c -p s f , zt] 
IF (a s > ac + ps) THEN 

2 l»gh 

V in =V in + \dzA in = Vin + AVs VIILA.4(9) 

where A V3 is given in the next section below. 
ELSE 

High 

K=V in + \dzA in = Vin + AV 4 VIILA.4(10) 



where A V4 is given in the next section below. 
END IF 

IF (zhigh = zt) keepGoing = FALSE 



15 



END IF 



END IF 



END IF 



! End Subcase 4.4 



! End Subcases 



! End hemisphere IF's 



Zlow — Zhigh 

END DO 



! Prepare for next z interval 
! End DO WHILE loop 
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2 hlgh 

Evaluation of the Integral \dzA in 



We now address the integral 



High 



dzA„ 



VIILA.4(11) 



where Am is given by Eqs. VIILA.4(2d, 3 ? and 4) as 



A. - a„ cos 



p*+a 2 c -a 2 s +C 



■2\ 



-1 Ps 



+ (^ 2 -^ 2 )cos 



-l 



+(P S +a e f -a)\a]-{p s -a c f 



VIII.A.4(12) 



with 



£ = z - zs . VIILA.4(13) 
This integral will involve elliptic integrals of the first, second, and third kinds 

AO which for present purposes are defined as in Abramowitz and Stegun. The 

jf incomplete elliptic integral of the first kind is given by 



F{ep,k)=\dG . 1 
o -vl — k si 



sin 2 0 



VIII.A.4(14) 



The complete elliptic integral of the first kind is denoted by K(k) and is equal to 
F{n/2,Tz). Similarly, the incomplete elliptic integral of the second kind is given by 



15 E{<p, k) = \de^\-k 2 sm 2 0 



VTII.A.4(15) 
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with the complete elliptic integral of the second kind denoted by E(k) and equal to 
E{n/2,k). Finally, the incomplete elliptic integral of the third kind is given by 



with the complete elliptic integral of the third kind denoted by Tl(n,k) and equal to 
5 Yl(n/2,n,k). Various efficient routines for the evaluation of these special functions 
can be found in the literature. 

Focusing on the integral of each of the three terms in Eq. VTII.A.4(12), it is 
possible to perform an integration by parts of the first term of Eq. VIILA.4(12) to 
dispatch the inverse cosine term. 



U(q> 9 n 9 k) = \d0 



(1-Hsin 2 ff)<Jl -k 2 sin 2 0 " 



VIILA.4(16) 



0 



WO a\ ]d£ 



cos 



[2 .2 2 . 

A +a c ~ a s 



) 



= a 2 c £cos 




+ 2a 2 c \d£ 



VIILA.4(17) 



where 




VIILA.4(18) 



The right-most integral of Eq. VIILA.4(17) is an elliptic integral. 



15 



An integration by parts also may be used for beginning the integration of the 



second term of Eq. VIII.A.4(12). 
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^ „2 „1 , 2 

Ps ~ a c + G s ~Q 



1 2 



= 4 a, — 



■1\ 



COS 



-1 



p)-al+a\-e 



•2 A 



3J 



cos 



+ \a 2 s {<-p 2 s))d^- 



1 



■2 \ nl/2 



VIII.A.4(19) 

where the integrals after the second equal sign in Eq. VIILA.4(19) come from 
1 partial fraction separation of the integrand after the first equal sign. Each of the 
i 5 integrals involves one or more elliptic integrals. 

;j The third term of Eq. VIII.A.4(12) is already in a form amenable to solution 

by elliptic integrals: 



VIII.A.4(20) 



4 Combining these three expressions in Eqs. VIILA.4(17, 19, and 20), yields the 
10 general expression for the AV used in Eqs. VIII.A.4(7 to 10): 



AV = a 2 c gcos 
2 



2 Ps a c 



1 



+ {[a]--?\cos 



p]-a]+a]-e 



+ - A 2 )A> +1 'A ~*\ + \p]\h +\h-\a]{al -p 2 s )P--E 



2 



VIII.A.4(21) 
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where 

r -\dC— VIII.A.4(22) 

i R 

L - he-?- VIILA.4(23) 

-^2 \ a h D l/2 ' 
S R 

L J\dr£- VIII.A.4(24) 

<L>4 — J u h D l/2 > 

i R 

P - \ dC I VIII.A.4(25) 

and 

E = )d£R V2 VIII.A.4(26) 

with R m given by Eq. VIII.A.4(18). The explicit solutions to the L, P, and E 
integrals of Eqs. VIII.A.4(22-26) depend on the values of £r, £?, a $ , and (a c - ps). 
These expressions are used in the preferred implementation. 

Explicit Solution for Interior Volume Vin When a s >ac+ps 
In the case where a s 2 > (a c + p s ) 2 , the preferred implementation follows the C, 
integration up the z axis, the southern hemisphere's sphere circle eventually 
circumscribes the cylinder circle unless it hits the top of the cylinder first. The 
limits of integration are therefore given by 

<Zi = MAX[B,-(zt-Zs)] and £2 = MIN[A, -(zb - z s )] VIII.A.4(27) 
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where A and B are given below. If a s 2 > (a c + fs) 2 , then it follows that a s 2 > (a c - fk) 2 . 
Therefore, we define 



-l-A' -<«.-«*■ Vin.A.4(28a) 

and 



B = ^-{a e+ p M f VIII.A.4(28b) 
and note that A > B > 0. Then 



R in = Jtf-tfXf-B 2 ) . VIII.A.4(29) 
The integrals for this form of R m are readily found in references such as 
Gradshteyn and Ryzhik. Define 



-i M 2 -^ 2 VIII.A.4(30a) 



^ =sin "Wl^F' VIII.A.4(30b) 
and 



^_ 4A 2 -B 2 ^ VIII.A.4(30c) 

A 

Then, recalling the definitions for the elliptic integrals of Eqs. VIII.A.4(14 to 16), 

the following expressions for the integrals are obtained: 

Lo = A-i [F(2i,g) - F(A 2 ,q)] , & < A VIII.A.4(31) 

Note that F(fa,q) = 0 if £ = A. 
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L 2 = A \E(M,q) - E(A*qr)] , & < A VIII.A.4(32) 

Again, EU 2 ,g)] = 0 if & = A. 

= .4/3 [2(A* + -B 2 ) E(Ai,g) - B» F(fc,g)] + £/3 [(A* - - B 2 )] 112 - 

A/3 [2{A2+m E(X 2 ,q)-B 2 F(/b,<?)] + &3 [(A^^i^-B 2 )] 1 ' 2 , 

VIII.A.4(33) 

Note that the entire second line is zero if £ = A. 



P = 



1 



A 2 -B 2 A ( A2 d2 A 



A -a, ) 



-n 



4 



VIII.A.4(34) 



The second elliptic integral withA;? as its first argument is zero if £2 - A. 
3 Lastly, 

Jo E = A/3 [(A 2 + B 2 ) E(Xi,q) - 2B 2 F(Ai,q)) - &I3 [(A 2 - - B 2 )] 112 ~ 

3 A/3 [(A 2 +B 2 ) E(vb,q)-2B» F(A 2 ,q)] - &I3 [(A*-£*)(0»-B 8 )P, ^<A 

VIII.A.4(35) 

'\ Again, the entire second line of this expression is zero if - A. 

% Explicit Solution for Interior Volume Vin when a s <a c +ps 

15 In the case where a s 2 < (ac + ps) 2 , A 2 is still defined as it was before, A 2 = a s 2 

(oc - ps) 2 , but B 2 becomes the negative of its previous definition: 
B» =(* + /*)*-*». VIII.A.4(36) 
Although A 2 is always greater than zero, we simply have that B 2 > 0. 
The radical is defined differently as 



20 R V2 = -yI(A 2 -£ 2 )(£ 2 +B 2 ) 



VIILA.4(37) 
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and that changes the elliptic integral solutions of Lo, L2, L4, P, and E. Defining 



y x = sin" 



y 2 = sm 



\Q U 2 +B 2 



AiB 2 +£ 



VIII.A.4(38a) 



VIII.A.4(38b) 



r = 



^A 2 +B 2 



VIILA.4(38c) 



the following expressions are obtained for the integrals. Note that £1 is inside 
absolute value signs so that only its magnitude is used. The upper limit £2, 
however, may be positive, negative, or zero. A negative £2 results in a change in 
sign for all terms in the integrals involving £2. All the integrands are even in ^so 
that all the integrals are odd. The elliptic integrals are also odd in their first 
argument. The expressions for the desired integrals are now 



1 



[F(y x ,r) + F(y 2 ,r)-\ 



VIII.A.4(39) 



L 2 =ylA 2 +B 2 E(r u r)- 



1>2 A 2 - £ 2 



4a 2 + b' 



^A 2 +B 2 E(y 2 ,r)- 



B 2 



4a 2 + b' 1 
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VIII.A.4(40) 
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L 4 = ■ 1 [(2B 2 -A 2 )B 2 F(y„r)-2(B 4 -A 4 )E( ri ,r)]- 
3^A 2 +B 2 

l f(2A 2 -B 2 + ^E^ 

3 * B +Q VIII.A.4(41) 

+ . 1 [(2B 2 - A 2 )B 2 F(y 2 ,r)-2(B 4 -A 4 )E(y 2 , r )] - 
3<JA 2 +B 2 



P = 1 ■ [B 2 U(r v n, r) + a 2 s F(y v r) + B 2 U{y 2 ,n,r) + a 2 F(y 2 ,r)] 

a 2 (a 2 s +B 2 )jA 2 +B 2 



VIII.A.4(42a) 



5 where 



A 2 (a 2 +B 2 ) VIII.A.4(42b) 



n ~a 2 (A 2 +B 2 ) 



and 



EAiA Y 7¥[B 2 F(y v ry(B 2 - A 2 )E( ri ,r)] + l -f(g +2B 2 -A 2 )^-^ 



+ X -jA 1 ^lB 2 F{y 2 ,ry{B 2 -A 2 )E{y 2 ,r)\ + ^(g +™ 2 - A 2 )^^ 

VIII.A.4(43) 

10 which completes the solutions for the intersection of a sphere and cylinder. The 
explicit expressions for £/ and are given in Section VIII.A.4.a. 
These expressions are used in the preferred implementation to determine the 
intersection. 
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The task of determining the packing fraction also may be accomplished 
according to an aspect of the invention using a modified approach. 

Recall that the radial distance of the sphere center from the z axis by 
p s = [x s 2 + yfiu* . (2-3-2) 

5 

Also recall that there are four possible cases for this overlap problem. In Case 
1, the sphere is completely outside the cylinder so Vm = 0. In Case 2, the sphere is 
completely inside the cylinder so Vm = (4/3)7ta s 3 . In Case 3, the cylinder is 
completely inside the sphere so Vin = %a c \zt - zt). In Case 4, the sphere and cylinder 

Jo intersect each other and the expression for Vin is more complex as will be derived 

!i below. 

a Case 1: Sphere Completely Outside Cylinder 

The sphere is completely outside the cylinder 
I V in = 0 (2-3-5) 
115 only if one of the following conditions are met: 

If the sphere's north pole is lower than the cylinder bottom . . . 
z s + a s < z b (2.3-6a) 
Or if the sphere's south pole is above the cylinder top . . . 
zs-as>zt (2-3-6b) 
20 Or if the sphere is more than its radius away from the cylinder side . . . 
A - as > ac (2-3-6C) 
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Or if the sphere center is above the cylinder top and is more than its radius away 
from the upper cylinder edge. . . 

ips - a c ) 2 + (z s - zi) 2 > a s 2 and z s > zt and ps>a c (2.3-6d) 

Or if the sphere center is below the cylinder bottom and is more than its radius 

away from the lower cylinder edge . . . 

0* - a c ) 2 + (zb - zs) 2 > a s 2 and z s <zb and ps>a c . (2.3-6e) 

Case 2: Sphere Completely Inside Cylinder 

The sphere is completely inside the cylinder 
Vin = (4/3)7ta s 3 ( 2 - 3 " 3 ) 
only if all three of the following conditions are simultaneously met: 
If the sphere's north pole is no higher than the cylinder top . . . 
Zs + a s <zt (2.3-4a) 
And if the sphere's south pole is no lower than the cylinder bottom . . . 
zs - as > z b (2-3-4b) 
And if the sphere's equator nowhere exceeds the cylinder radius . . . 
ps + a s < a c . (2-3-4C) 

Case 3: Cylinder Completely Inside Sphere 

The cylinder is completely inside the sphere 
Vin = nac 2 (zt - Zb) ( 2 - 3 " 7 ) 
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if both the following two conditions are met: 
If the upper edge of the cylinder is inside the sphere . . . 

{ps + a c ) 2 + (zt - zb) 2 < as 2 (2.3-8a) 
And if the lower edge of the cylinder is inside the sphere . . . 

5 (ps + a c ) 2 + (z s - z b ) 2 < as 2 . (2.3-8b) 

Case 4: Cylinder and sphere intersect 

Again, there are several subsets of Case 4. To efficiently organize them, we 

define two circles: the sphere circle and the cylinder circle. The sphere circle is the 
10 perimeter of the sphere's cross section at a given z altitude. Its radius is 

r s = [as 2 -(z- zs) 2 ] 1/2 , zs - a s < z < z s + a s (2.3-9) 

and its radial position ps is given by Eq. (2). Likewise the cylinder circle is the 

perimeter of the cylinder's cross section which has constant radius a c at any z 

altitude between zb and zt and is centered on the z axis. 
15 To describe the intersection of sphere and cylinder, we define four subcases 

for he relationship of the sphere circle and the cylinder circle at a particular 

altitude z: 

Subcase 4.1: The sphere circle is completely enclosed in the cylinder circle. 
Subcase 4.2: The two circles do not overlap or one or both are nonexistent. 
20 Subcase 4.3: The cylinder circle is completely enclosed in the sphere circle. 
Subcase 4.4: The two circles intersect. 
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To find Vin, our task is reduced to integrating in z the overlap areas of the two 
circles Ain from the altitude where one subcase holds to the altitude where another 
subcase holds. 

5 For Subcase 4. 1, An = %a s 2 . (2.3-10a) 
For Subcase 4.2, Ain = 0. (2.3-10b) 
For Subcase 4.3, An = na<?. (2.3-10c) 
For Subcase 4.4, Am = a c 2 {a - cosa since) + r s 2 (J3 - cos/? sinfi) (2.3-10d) 
where r s is given by Eq. (2-3-9) and 

10 cosa ^*"'-'- ( 2 - 3 - lla ) 



^ ^ 2 -K-A) 2 ][(^+A) 2 -^ 2 ] (2 . 3 . nb) 

2 P, a c 

cos/3= P ° +r *~ a < (2.3-128) 
IPS, 



rin^- V^-^-A ^+A) 2 -^] (2 .3-12b) 

It is now necessary to determine a series of altitudes at which the cylinder 
15 circle and the sphere circle change from one subcase to another in Case 4. 

A sphere and a cylinder side can intersect each other at either 0, 2, or 4 
altitudes. Pairs of altitudes may be degenerate (two solutions having the same 
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altitude) if the sphere just nicks the cylinder wall. The solutions for these altitudes 
come from solving the following equation: 

Mz)] 2 =[ac±/*] 2 (2.3-13a) 
or 

5 as 2 -(z-zs) 2 =[ac±ps] 2 (2.3-13b) 
which, using Eq. (2-3-9), has the four solutions 

z = z s ±^-(a c ±p s f . (2-3-14) 

The + sign in the first ± of Eq. (2-3-14) denotes the upper or northern hemisphere 

intersection altitude; the - sign denotes the lower or southern hemisphere 

10 intersection altitude. The + sign in the second + of Eq. (2-3-14) denotes intersection 
with the far cylinder wall, on the opposite side of the z axis from the sphere's side; 
the - sign denotes intersection with the near cylinder wall, on the same side of the z 
axis as the sphere. 

Armed with these four solutions, we will now consider each of the four 

15 subcases of Case 4, but we will consider separately whether the sphere circle lies in 
the southern hemisphere of the sphere or whether it is the equatorial plane or 
northern hemisphere. We consider the four subcases separately for the two 
hemispheres making eight situations in all. We consider the hemispheres 
separately because if the sphere circle is in the southern hemisphere, it is 

20 increasing in size as z increases. Hence the preferred implementation can 
determine the altitude at which subcases change. If the sphere circle is its 
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equatorial plane or in its the northern hemisphere, it is decreasing in size as z 
decreases and the altitude at which subcases will change again can be determined. 

We now outline the overall logic used in the preferred implementation to find 
5 the altitudes at which the subcases for Case 4 change. Let ziow be the lower altitude 
of the subcase region currently being examined and zu 8 h be its upper altitude. 

The absolute bounds below or above which there can never be any overlap 
between sphere and cylinder is given by 

z 0 =MAX[z b ,z s -a s ] (2-3-15) 

10 for the lower bound and 

z 5 = MJN[z t ,z s + a,] (2.3-16) 
for the upper bound. The altitudes z\ through za correspond to the different 
intersection altitudes (if they exist) of the sphere's polar plane that is perpendicular 
to the cylinder axis. A sphere's polar plane is the plane which intersects the sphere 

15 in a great circle perpendicular to the equatorial plane. It can be oriented at any 
azimuthal angle. The polar plane of interest to us is the one passing through the 

cylinder axis. The intersection altitudes are illustrated in Fig. . 

In ascending order (again, if they exist) the intersection altitudes are: 
Southern hemisphere intersection with cylinder wall closest to sphere center 

20 z x =z s -^ s -{a c -p s f (2.3-17) 

Southern hemisphere intersection with cylinder wall furthest from sphere center 
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Northern hemisphere intersection with cylinder wall furthest from sphere center 
z,=z s+ ^-(a c+Ps f (2-3-19) 
Northern hemisphere intersection with cylinder wall closest to sphere center 

5 * 4 =z,+V«. 2 -te-A) 2 . < 2 - 3 - 20) 

In all four of these solutions, if the argument of the square root is negative, the 
square root is replaced by zero. If the sphere intersects both sides of the cylinder, 
there are four real solutions; if only one side, there are two real solutions; if not at 
all, there are no real solutions of Eqs. 17 through 20. The figure below shows the 

10 two real solutions for a sphere that intersects only one side of the cylinder. 

The method used in the preferred implementation for finding the volume Vin 
of a sphere that lies inside a cylinder will now be discussed. The method is divided 
into two DO WHILE loops, one for the southern hemisphere and one for the 
northern hemisphere (which includes the sphere's equatorial plane). This is 

15 convenient because for the southern hemisphere the sphere circle is always 
increasing with increasing z while for the northern hemisphere it is always 
decreasing. The method as specifically implemented here has a number of 
comments so that the reader can follow the reasoning guiding it. 



20 Vm = 0 ! Vin is the volume of this sphere lying inside the 
! cylinder. Initialize it to zero, 
zo = MAX[zb, z s - a s ] ! z 0 and zs are the absolute lower and upper 
zs = MIN[zt, z s + a s ]! limits of possible sphere-cylinder overlap. 
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IF (z s > Zb) THEN ! The southern hemisphere can contribute to Vi n if and 
! only if this statement is true. Otherwise it lies 
5 ! beneath the cylinder bottom and cannot contribute. 

! (It might not contribute even when the above 
! statement is true but that is checked below.) 

Zlow = Zo ! Initialize zi ow to smallest possible value. 
10 r s (ziow)=[a s 2 -(ziow-z s ) 2 ] 1/2 ! Find initial sphere circle 

! radius. It is 0 if square 
! root argument < 0. 



15 ! Find the initial southern hemisphere subcase. 
IF (p s + r s (ziow) < a c ) THEN 

subcase = 1 ! Sphere circle inside cylinder circle. 

ELSE IF (p s - r s (ziow) > a c ) THEN 

subcase = 2 ! Sphere circle outside cylinder circle. 
20 ELSE IF (r s (ziow) - ps > a c ) THEN 

subcase = 3 ! Cylinder circle inside sphere circle. 

ELSE 

subcase = 4 ! Sphere circle intersects cylinder circle. 
END IF 

25 

! We exit the DO WHILE loop when the subcase is no longer a positive 
! integer. Remember that since we are in the southern hemisphere, the 
! sphere circles are expanding with increasing z. 
30 DO WHILE (subcase > 0) ! Southern hemisphere loop. 

IF (subcase = 1) THEN 

! The sphere circle lies inside the cylinder circle. 

! It expands with increasing z until one of three 
35 ! things happen: 

! (1) It hits the top of the cylinder. 

! (2) It intersects the cylinder circle. 

! (3) It reaches its equatorial plane. 

! Numbers (2) and (3) can be combined if we agree to 
40 ! define zi with the square root replaced by zero if 

! its argument is negative. 



w = a s 2 - (a c - ps) 2 
IF (w > 0) THEN 
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Zl = Z s - W 1/2 

ELSE 

Zl = Z s 

END IF 

Zhigh = MIN[zt, zi] 

Vin = Vin + Vi(ziow,Zhi g h) ! We'll derive an explicit 
! expression for this Vi later. 



! Now we must find the next subcase. 
10 IF (zhigh < MIN[z s ,z t ]) THEN ! There is another 

! subcase to be 
! addressed in the 
! southern hemisphere. 
IF (p s > 0) THEN ! Sphere is not concentric with 
15 subcase = 4 ! the cylinder so it will 

! expand until it intersects. 
ELSE ! Here, the sphere is concen- 

subcase = 3 ! trie so it will expand 
END IF ! until it perfectly circum- 

20 ! scribes the cylinder 

! circle. 

ELSE ! You reached the end of southern 

subcase = 0 ! hemisphere contribution. End the 
END IF ! S. hemisphere's DO WHILE loop. 

25 

ELSE IF (subcase = 2) THEN 

! The sphere circle lies outside the cylinder circle 
! until it expands to the point where it intersects 
30 ! the cylinder. This intersection must happen or we 

! wouldn't be here in Case 4. The intersection 
! happens at zi. There is no contribution to Vi n from 
! this region. We only need to find zhigh = zi. 

35 w = a s 2 - (a c ■ ps) 2 

IF (w > 0) THEN 

Zhigh - Z S - W 1/2 

subcase = 4 ! Not done yet with S. hemisphere. 

ELSE 

40 Zhigh = z s 

subcase = 0 ! Reached equatorial plane. Now 
END IF ! we're done with S. hemisphere. 
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ELSE IF (subcase = 3) THEN 

! The cylinder circle lies inside the sphere circle. 
! Since the sphere circle is expanding, it will 
! remain inside until it reaches the first of either 
! the top of the cylinder or its equatorial plane. 

zhigh = MIN[z t , z s ] 

Vin = Vin + V 2 (ziow,Zhi g h) ! We'll derive an explicit 
! expression for this V2 later. 

subcase = 0 

ELSE IF (subcase = 4) THEN 

! The only possibility left is that sphere and 

! cylinder circles initially intersect. This region 

! extends upward until one of three things happen: 

! (1) We reach the top of the cylinder. 

! (2) We reach the sphere's equatorial plane. 

! (3) We reach a point at which the sphere circle 
! completely circumscribes the cylinder circle. 

! We cannot reach the top of the sphere in this case 

! because we are only dealing with the southern 

! hemisphere. (Remember Z2 = z s if the argument of 

! the square root is negative.) 

w = a s 2 - (a c + ps) 2 
IF (w > 0) THEN 

Z2 = Z s - W 1/2 

subcase = 3 

ELSE 

Z2 = Zs 

subcase = 0 
END IF 

Zhigh = MIN[z t , z s? Z2] 

Vin = Vin + V 3 (ziow>zhi g h) ! Well derive a numerical 
! integration for this V3 later. 

END IF 

! If you have reached either the cylinder top or the 

! sphere's equatorial plane, you are done with the southern 

! hemisphere's contribution to Vin. Go to the northern 

! hemisphere next. (We use > signs below in case of round- 

! off error. With infinite precision, Zhigh could never be 

! greater than zt or z s .) 
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IF (zhigh > z t OR Zhigh > z s ) THEN ! Redundant check saying 

subcase = 0 ! we're done with the southern hemisphere. 

ELSE 

ziow = zhigh ! Cycle back through the DO WHILE loop 
5 END IF ! for the next region in the southern 

! hemisphere contributing to Vin. 
END DO ! End of southern hemisphere DO WHILE loop. 

END IF ! End of southern hemisphere overall IF statement. 
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IF (z s < z t ) THEN ! The northern hemisphere can contribute to Vin if and 
! only if this statement is true. Otherwise it lies 
15 ! above the cylinder top and cannot contribute. (It 

! might not contribute even when the above statement 
! is true.) 

ziow = MAX[zb, z s ] 

20 r s (ziow)=[a s 2 -(ziow-z 8 ) 2 ] 1/2 ! Find initial sphere circle radius in 

! N. hemisphere. It is 0 if square root 
! argument < 0. 

! Find the initial northern hemisphere subcase. 
25 IF (p s + r s (ziow) < a c ) THEN 

subcase = 1 ! Sphere circle inside cylinder circle. 

ELSE IF (p s - r s (ziow) > ac) THEN 

subcase = 0 ! Sphere circle outside cylinder circle. 

! Although this is actually subcase 2, since 
30 ! N. hemisphere sphere circle can only contract 

! it will never intersect the cylinder 
! circle if it does not initially intersect it. 
! Therefore, we set the subcase to zero so the 
! N. hemisphere's DO WHILE loop will end. 
35 ELSE IF (p s + ac < r 8 (ziow)) THEN 

subcase = 3 ! Cylinder circle inside sphere circle. 

ELSE 

subcase = 4 ! Sphere circle intersects cylinder circle. 
END IF 

40 

DO WHILE (subcase > 0) ! Begin northern hemisphere loop. 



IF (subcase = 1) THEN 

140 



! The sphere circle lies inside the cylinder circle. 
! Since the sphere circle shrinks as z increases, it 
! cannot ever intersect the cylinder circle. It can 
! only end at its north pole or at the top of the 
! cylinder. That ending position is zs. 

Zhigh — Z5 



V in = Vin + Vi(ziow,zhi g h) ! We'll derive an explicit 



subcase = 0 



ELSE IF (subcase = 2) THEN 

! The sphere circle is inside the cylinder circle. 
! In the northern hemisphere this means they will 
! never intersect more so there will be no more 
! contribution to Vin. 

subcase = 0 



ELSE IF (subcase = 3) THEN 

! The cylinder circle lies inside the sphere circle. 
! The sphere circle will contract until we hit top 
! of the cylinder or until we reach Z3 which is the 

! point at which the two circles will intersect. 

w = a s 2 - (a c + p s ) 2 ! w is nonnegative or we 

! couldn't be here. 

Z3 = Z s + W 1/2 

Zhigh = MIN[z t , zs] 

Vm = Vm + V 2 (ziow ? Zhi g h) ! Well derive an explicit 
! expression for V2 later. 

IF (zhigh < z t ) THEN 

IF (p s > 0) THEN ! Sphere is not concentric 
subcase = 4 ! with cylinder so it will 
! contract until it inter- 
! sects the cylinder circle. 



! expression for Vi later. 



ELSE 



! The sphere is concentric so 



subcase 



= 1 ! it will contract until it 



END IF 



! suddenly lies completely 



ELSE 



subcase = 0 



! inside the cylinder circle. 
! If we're here, we reached 
! top of the cylinder; we're 
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END IF 



! done with N. hemisphere. 



ELSE IF (subcase = 4) THEN 

! The only possibility left is that sphere and 
! cylinder circles initially intersect. This region 
! extends upward until one of two things happen: 
! (1) We reach the top of the cylinder before 
! reaching the north pole. 
! (2) The intersection stops with the sphere circle 
! now either completely inside or completely 
! outside the cylinder circle. This happens at 

! Z4. 



w = a s 2 - (ac - ps) 2 ! w is never negative or we 

! wouldn't be here. 

Z4 = Z s + W 1/2 

Zhigh = MIN[z t , z 4 ] 

Vin = V in + V 3 (ziow,Zhigh) ! We'll derive a numerical 

! solution for V3 later. 

IF (p s < ac) THEN 

subcase = 1 ! Next, sphere circle will lie 
! inside cylinder circle. 

ELSE 

subcase = 0 ! Next, sphere circle will lie 
! outside cylinder circle. 

END IF 



END IF 

! (Note the second logical expression in the IF statement below is z 5 , 
! not z s .) 

IF (zhigh > z t OR zhigh > z 5 ) THEN ! Redundant check saying 
subcase = 0 ! we're done with northern hemisphere. 

ELSE 

ziow = Zhigh ! Cycle back through the DO WHILE loop 
END IF ! for the next region in the northern 

! hemisphere's contribution to Vin- 

END DO ! End of northern hemisphere DO WHILE loop. 

END IF ! End of northern hemisphere overall IF statement. 
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Expressions for the three volume expressions Vi(ziow, zhigh), i = 1, 2, or 3 will 
now be addressed. Remember in the following that a c is the cylinder radius, a s is 

the sphere radius, (x s , y s , z s ) is the position of the sphere's center, and p s = ^x] 

is the radial distance of the sphere center from the cylinder axis. 

5 

Derivation of Vi. 

The volume V\ is an integral of full sphere circles from ziow to zhigk because the 
sphere circles lie entirely within the cylinder circle. 

2 high 2 high 

V x =7t \dz[r s {z)f=n ]dz[a'-(z-z s ) 2 ] (2.3-21) 

z }cw 2 law 

10 or 

V x = j(z higk -zjftf -z 2 high -zl+z^-zJ-Sz^-zj] (2.3-22) 

where ziow and zhigh depend on where we are in the above method. Their values are 
given in the implementation shown above. 



15 Derivation of V2. 

The volume V2 is the integral of the cylinder circle from ziow to zugh because 
the cylinder circle lies entirely within the sphere circle. It is given simply as 

V 2 = 7ra 2 c {z high -z low ). (2.3-23) 
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Derivation of Vs. 

The volume V3 is the integral of intersecting sphere and cylinder circles from 
ziow to zhigh. The limits are always given in the method as described above. The 
intersect area is given by Eqs. 2-3-10, 11, and 12 which we must integrate in z. 

5 V % = jdz{al[a-cosasma] + [a^-(z-z s ) 2 ][jB~cosj3sin^]} (2.3-24) 

z low 

where 



cos a = 



P\ , +fl e -fl;+(z-z,) 
2PA 



(2.3-25) 



sma = 



_ VK 2 - («, - Ps f-{z~z s ) 2 ][(a c + A) 2 ~ "s + ~ *,) 2 ] 



(2.3-26) 



a = arccos 



(2.3-27) 



10 and 



2p sA /« 2 -( z - z s ) 2 
an I = ^ " ( *' ~ A)2 " (Z ~ Z ' )2][(fl ' + Ps)2+ *' ~ (z ~ ^ 



2p s ^a 2 s -(z-z s f 



(2.3-28) 



(2.3-29) 



(3 -- arccos 



A 2 -« 2+g2 -( z -^) : 



(2.3-30) 



The integral of Eq. 2-3-24 could be done analytically in terms of complete and 
15 incomplete elliptic integrals of the first, second, and third kinds. However, this 
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solution is so cumbersome that it is more efficient to perform it numerically, 
especially because the integrand is quite smooth and well-behaved so that a 
Romberg integration is very efficient and contains an accurate error assessment. 



5 Calculating the Total Volume Fraction for a Pack 

The preferred implementation also can determine the volume fraction § of the 
pack. Only the innermost cylinder contains particles that look like a full three- 
dimensional pack. The outer cylinders do not contain the particles the inner 
cylinders contain. However, it would be unfortunate to throw away all the pack 

10 information the outer cylinders contain. They do perturb the pack of the innermost 
cylinder correctly and their correct geometry for the particles they contain should 
somehow be used in estimating the volume fraction of the pack. 

The preferred implementation determines the volume fraction 0in the 
following way. The j-th cylinder can contain all particle modes whose cylinder radii 

15 are less than or equal to its cylinder radius a c (j). As long as the water level option 
of pack building is operating, the j-th mode particles are distributed approximately 
as they would be in a truly three-dimensional pack. The preferred implementation 
therefore calculates the volume fraction $ of Mode j particles within the j-th 
cylinder and calls it the j-th partial volume fraction. It does likewise for all the 

20 other modes, e.g., J of them. Then the total volume fraction is simply the sum of all 
the partial volume fractions: 
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* = 2>,. VIILA.5(1) 

If desired, the user can specify that only the volume fraction of the innermost 
cylinder, where the pack is three dimensional, will be calculated. 

To calculate the partial volume fraction $ in the j-th cylinder, the preferred 
5 implementation sums the interior volume of each Mode j sphere inside the j-th 
cylinder. Let rtj be the number of Mode j spheres in the j-th cylinder. Let Vj(k) be 
the volume inside the cylinder of the fe-th sphere of Mode j. Then 

V<o t U) = £vj(k) • VIILA.5(2) 

Dividing this total interior partial volume by the volume of the j-th cylinder VcyiQ) 
10 gives the partial volume fraction $ for the pack: 

4i = Vtot(j)/V C yi(j) VIILA.5(3) 
where 

VcyiQ) = n[a c (j)nztop(j) - z bo t(j)] . VIII A. 5(4) 

The user has the option of making calculating volume fractions in cylinders 
15 that are a little smaller than the actual cylinders used in making the pack so that 
possible wall effects can be eliminated from the calculation. Remember, however, 
that the cylinder walls are barriers to the particle centers, not their leading edges 
so that disturbances to the pack from wall effects normally will be small. 

20 Finding the Radial Distribution Functions 
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The radial distribution functions, denoted by gij(r)dr> are probability 
distribution functions denoting the mean number of Mode or Submode j particles 
whose centers are between r and r + dr from the centers of Mode or Submode i 
particles throughout the pack. They can be important in establishing many 
5 properties of random packs. Because the packs here are assumed random, only the 
scalar radial distance r is used rather than vector r in this implementation. 

Reduced dimension packs do not extend equally in all directions. They are a 
series of concentric cylinders. Therefore, for a given Mode or Submode i particle, 
the preferred implementation calculates the volume of the spherical shell segment 

10 centered on Particle i and lying between r and r + dr, that lies inside the j-th 

cylinder. Then it counts all Mode j particles lying in that shell segment. Dividing 
the number of Mode j particles by the volume of the shell segment gives an 
approximation to gij(r)dr at that position r. 

Finding the volume of this spherical shell segment is not a trivial procedure. 

15 Fortunately, a means to calculate the total sphere volume lying inside any cylinder 
is described above. Therefore, the preferred implementation can calculate the total 
volume of an imaginary sphere of radius r + dr, centered at Particle i, that lies 
inside the j-th cylinder, and subtract from it the total volume of the concentric 
sphere of radius r to obtain the volume of the spherical shell segment between r and 

20 r + dr lying in the j-th cylinder. 
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The user specifies the maximum radial distance r m ax(ij) for which gij is 
desired and the number of bins mij into which it divided rmax(ij). Then the bin 
widths have finite width 

Anj = rmax(ij)lmij VIILB(l) 
5 rather than infinitesimal width dr. Then the gij are determined as follows. 

First, note that the gij(r) are zero for r<ai + aj. Then, consider a particular 
Mode i particle with index number ih. Denote the number of Mode j particles in the 
first bin from Particle ih particle by Nj{ih\ at + aj). The first bin lies between (a; + aj) 
and (at + aj) + An/. In general, the number of Mode j particles lying in the k-th bin 
10 from Particle ih is denoted by Nj(h; rk) where rk = ai + ay + (k-l)Anj, k = 1, 2, ... , mij. 

Now consider the volume of the k-th bin from Particle ih that lies inside the jf- 
th cylinder. This can be obtained from the volumes Vin of previous sections where 
Vin = Vin(x s , ys, Zs, a s ; a c , zbot, ztop). The argument list contains the position of the 
sphere center (xs, ys, z s ), the sphere radius a$, the cylinder radius a c , and the 
15 cylinder's bottom and top altitudes (zbot, ztop). The volume for the k-th bin about 
Particle ih is then given by 

Vin(ih,j; k) — Vin(xih, yih, Zih, Vk+1\ acj, Zbotj, Ztopj) — Vin(xih, Jih, Zih, Tk\ acj, ZbotJ, Ztopj) 

VIII.B(2) 

and the number density of Mode j particles in this volume is 
20 o{ih,j;k) = Nj(i h ; r k )/Vin(ih,j; k) . VIII.B(3) 
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Finally, to obtain the gij(r), the above expression can be averaged over all 
Mode i particles: 

vm.B(4) 

if r > ai + aj and is zero otherwise. The term k is the integer lying in the range 

5 

(r-ai- aj)/Anj < < 1 + (r-ai- aj)/Anj . VIILB(5) 
If r is so large or so small for a particular particle ih that the spherical shell 
containing it misses the j-th cylinder altogether, then that term in the sum of Eq. 
VTII.B(4) is omitted altogether. 
10 The normalized radial distribution function, denoted Gij(r), is given by 



It is the mean number density nj of Type j particles at a distance r from a Type i 
15 particle. Thus normalized, the function Gy(r) approaches the constant nj as r grows 
large rather than increase quadratically as gij(r) does, nj represents the number of 
j particles per unit volume, e.g., concentration, output from the preferred 
implementation. 



Gij(r)=gij(r)/(4nr2) 



VIILB(6) 



or in our finite bin width treatment, it is given by 




VIII.B(7) 



4 
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Finding the Coordination Numbers 

The coordination number Qij is the mean number of Mode; particles touching 
a Mode i particle. Again, because a reduced dimension pack is ordered in concentric 
cylinders rather than being isotropic in all directions, the preferred implementation 
finds how much of the surface of the i-th particle is available for Mode ; particles to 
touch. That is equivalent to asking how much of the imaginary spherical surface 
centered on the Mode i particle and having radius at + aj lies inside they-th 
cylinder. We will denote that interior portion of the surface by Sij. Then if a search 
shows there are hij Mode ; particles touching a particular Mode i particle identified 
by index number ih, then our best estimate of the total number of touching Mode; 
spheres, if the entire surface were accessible, would be 

m = hij Mm + ai)ySij . VHLC(l) 

If the entire surface is accessible, then Sij = 4n(cn + aj) 2 and rjij = hij. Because 
of round-off error, a certain center-to-center error tolerance s is required to 
determine whether two spheres are indeed touching. The code considers them 
touching if their centers are within the distance (a; + a;)(l + s) from each other. 
Because the preferred implementation is in double precision, the default for s is 
10~ 12 , but the user can specify any other desired 8. Normally it should be a few 
orders of magnitude greater than machine precision. 

The spherical surface Sij lying inside a cylinder could be obtained from our 
earlier work on obtaining the spherical volume inside a cylinder by taking a 
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derivative of the volume expressions with respect to the sphere radius. However, 
those expressions are complicated and have substantial bookkeeping associated 
with them. Because in this implementation there is an error tolerance, we may 
take a numerical derivative as follows: 
5 Sij = e~ l [Vin(xih, yih, Zih, (at + aj)(l + s); a c (j), zbotj, ztopj) - 

Vin(xih, yih, zih, (ai + aj)\ a c (j), zbotj, ztopj)] . VIILC(2) 
The s used in this expression should be kept very small even if the user chooses to 
be rather liberal in assigning a rather large error tolerance. One might, for 
example, assign this s to be either the user-defined (or default) error tolerance or 
10 1CT 6 , whichever is smaller. 



Finding Various Pack Statistics 

There are a number of vital statistics of the pack that will generally be of 
interest. If the user so specifies, then the following data is written to a user-named 
15 file. 

1. A recapitulation of the particle histogram. 

2. A recapitulation of the values for control flags and parameters used in 
building the pack. 

3. The total number of particles placed in the pack and the number placed 
20 from each mode. 
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4. The minimum and maximum radial and axial positions for the centers of 
all particles in the pack and the same for each mode. 

5. Radial and axial particle density functions for each mode specified 
between the minimum and maximum positions specified above. 

6. The number of particles placed after tunneling (if tunneling is permitted). 

7. The number of particles hitting one of the catch nets. 

8. The number of particles hitting one of the water levels. 

9. The number of preplaced particles, if any. 



10 Additional advantages and modifications will readily occur to those skilled in 

the art. Therefore, the invention in its broader aspects is not limited to the specific 
details, representative devices and methods, and illustrative examples shown and 
described. Accordingly, departures may be made from such details without 
departing from the spirit or scope of the general inventive concept as defined by the 

15 appended claims and their equivalents. 
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